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.Abstract 

The  primary  topic  of  this  thesis  is  the  multigrid  solution  of  the  steady  state  2D  non¬ 
linear  conservation  law.  The  difficulties  encountered  in  the  numerical  solution  of  more 
complicated  problems  in  fluid  dynamics  (e.g.  the  steady  2D  Euler  system),  like  the  non- 
efliptidty  of  the  equations,  and  the  presence  of  discontinuities  in  the  solution,  can  be 
studied  in  this  model  scalar  case. 

The  work  deals  with  two  main  issues:  development  of  new  discretization  schemes  and 
adaptation  of  the  coarsening  techniques.' 

New  genuinely  2D  (based  on  a  “9-point  square”  stencil)  conservative  discretization 
schemes  are  developed.  These  schemes  provide  the  possiblity  to  separate  treatment  of  the 
streamwise  and  cross-stream  directions.  Due  to  this  separation,  the  artificial  viscosity 
can  be  added  in  the  streamwise  direction  only.  High  resolution  is  introduced  in  the 
cross-stream  direction.  Therefore,  the  resulting  schemes  have  good  stability  properties, 
are  second  order  accurate  and  provide  a  good  resolution  of  discontinuities:  representing 
them  in  the  numerical  solution  by  thin  oscillation-free  transition  layers. 

The  adaptation  of  the  coarsening  techniques  is  based  upon  a  more  precise  understand¬ 
ing  of  what  should  be  meant  by  discontinuity  location  in  a  shock-captured  solution.  We 
have  shown  that  a  such  solution  (provided  the  discretization  scheme  employed  is  conser¬ 
vative  and  seoond  order  accurate)  contains  information  about  the  discontinuity  location 
with  seoond  order  accuracy.  The  conventional  coarsening  techniques  (full-weighted  resid¬ 
ual  transfer  and  bilinear  correction  interpolation)  appear  to  provide  a  good  correction 
for  the  discontinuity  location  as  well  as  for  the  solution  in  smooth  regions. 

As  the  result  an  efficient  multigrid  solver  is  constructed  for  a  general  steady  state  2D 
conservation  law. 


Contents 


1  Introduction  1 

2  Conservation  law  and  its  discretization  4 

2.1  Boundary- value  problem:  solution  properties .  4 

2.1.1  Linear  convection  diffusion  equation .  4 

2.1.2  Nonlinear  equation .  5 

2.2  Discretization  of  conservation  law .  7 

2.3  Recovering  of  discontinuity  location .  7 

2.3.1  Integral  relation .  7 

2.3.2  Implementation .  10 

3  Multigrid  methods  12 

3.1  Elliptic  case . 12 

3.1.1  Relaxation .  12 

3.1.2  Coarse  grid  correction . 12 

3.1.3  Full  approximation  scheme .  13 

3.1.4  Multigrid  cycle .  14 

3.1.5  Full  multigrid  algorithm .  14 

3.2  Non-elliptic  case .  15 

4  Discretization:  Linear  case  16 

4.1  Some  existing  methods . 16 

4.1.1  Central  differencing .  16 

4.1.2  Upstream  differencing  .  17 

4.1.3  Upstream  second  order  scheme  .  17 

4.1.4  Remark  on  overshoots .  18 

4.1.5  Approach  of  Spekreijse .  18 

4.2  '2D  approach:  Homogeneous  equation .  19 


i 


ii  _ _ _  _ CONTENTS 


.  4.2.1  2D  scheme .  19 

4.2.2  First  order  N  scheme .  21 

4.2.3  S  scheme .  21 

4.2.4  Examples  of  limiters  .  .  .  . . 23 

4.2.5  Relaxation  . . .  .  .  . . . .  .  .  .  24 

4.3  2D  approach:  Inhomogeneous  equation . .  .  .  .  . .  25 

4.3.1  2D  scheme . '. .  25 

4.3.2  N  scheme .  26 

4.3.3  S  scheme .  27 

5  Discretization:  General  case  28 

5.1  Homogeneous  equation .  28 

5.1.1  Central  and  upstream  schemes .  28 

5.1.2  2D  scheme .  29 

5.1.3  First  order  N  scheme .  30 

5.1.4  Si  and  S2  schemes  .  . .  32 

5.2  Inhomogeneous  equation .  36 

5.2.1  2D  scheme .  36 

5.2.2  N  scheme .  38 

5.2.3  S  scheme .  39 

6  Numerical  experiments  41 

6.1  Linear  equation .  42 

6.1.1  Smooth  solution .  42 

6.1.2  Resolution  of  contact  discontinuities . 45 

6.2  Nonlinear  equation .  53 

6.2.1  Smooth  solution .  53 

6.2.2  Discontinuous  solution .  55 

6.2.3  Non-constant  solution  containing  discontinuities .  60 

6.3  Choice  of  discretization .  65 


CONTENTS 


tu 


7  Discussion  and  conclusions  66 

7.1  Summary .  66 

7.2  Remark  on  double  discretization . . .  .  66 

7.3  Extension  to  Euler  equations . - .  67 

7.4  Efficiency  comparison  .  .  .  .-.  .  . . ;  .  .  .  . . - .  67 

.  7.5  Some  properties  and  future  development  ............  ^  ....  .  69  . 


Chapter  1 
Introduction 


A  great  need  for  accurate  simulations  of  flows  with  discontinuities  exists  in  many  fields 
of  physics  and  engineering.  The  major  difficulties  in  this  area  are  the  non-elliptidty  of 
the  governing  equations  and  the  problem  of  representing  the  shocks  and  contact  discon¬ 
tinuities. 

Multigrid  methods  were  initially  developed  for  elliptic  problems  and  have  been  shown 
to  be  extremely  efficient  in  this  case.  However,  the  non-elliptic  case  still  requires  inves¬ 
tigation  in  order  to  achieve  the  same  efficiency. 

The  objective  of  this  work  is  to  develop  a  fast  multigrid  solver  for  the  scalar  steady 
state  2D  conservation  law.  Such  a  conservation  law  can  serve  as  a  good  model  problem 
for  more  general  non-elliptic  problems  (like  the  steady-state  2D  Euler  system).  This  is 
because  some  difficulties  encountered  in  their  numerical  solution,  such  as  the  conflict 
between  numerical  stability  and  higher-order  accuracy  or  the  presence  of  discontinuities 
in  the  solution,  can  be  studied  already  in  the  scalar  case.  Therefore,  this  work  can  be 
regarded  as  a  preparational  step  towards  developing  a  fast  solver  for  more  complicated 
non-elliptic  problems. 

This  solver  is  required  to  approximate  the  solution  such  a  conservation  law  fast, 
approaching  steady  state  directly  without  an  explicit  marching  in  time.  We  want  it  not 
only  to  produce  second  order  accurate  solution  in  smooth  regions,  but  also  to  deal  with 
discontinuities  in  some  efficient  way.  Therefore,  we  explore  the  problem  of  representing 
discontinuities  and  getting  them  to  converge  to  their  correct  position  by  a  multigrid 
solver  as  efficiently  as  obtaining  the  solution  in  smooth  regions  is  explored. 

One  approach  to  representing  discontinuities  by  a  finite  difference  method  is  to  use 
difference  equations  only  in  smooth  regions.  The  discontinuity  itself  is  followed  explicitly 
using  some  jump  conditions  supplemented  by  characteristic  data.  This  approach  is  called 
shock  fitting.  A  multigrid  fast  solver  can  be  used  both  for  obtaining  the  solution  in 
smooth  regions  and  for  following  the  path  of  a  discontinuity  in  an  efficient  way.  It  is  also 
possible  to  show  that,  once  we  obtain  higher  order  accuracy  in  the  smooth  regions,  we  can 
obtain  the  same  order  of  accuracy  in  the  discontinuity  location.  However,  for  complex 
flows  with  several  intersecting  shocks  such  procedures  become  difficult.  Furthermore, 
there  is  the  additional  difficulty  of  predicting  the  generation  of  shocks  that  are  not 
present  initially. 

Another  alternative  is  the  so-called  shock  capturing  method.  In  the  flow  of  a  real  fluid 
there  are  no  discontinuities.  There  are  instead  very  thin  regions  of  very  steep  gradients. 


Chapter  1  —  Introduction 


This  is  because  terms  representing  viscosity  and  heat  conduction  in  the  usual  equations 
of  fluid  dynamics  have  very  small  but  still  non-zero  coefficients.  However,  the  distance 
scale  on  which  the  resulting  transition  layers  axe  smooth  in  the  physical  solution  may  be 
smaller  than  any  reasonable  meshsize.  The  numerical  method  then  is  to  increase  the  size 
of  these  dissipative  terms  so  that  the  flow  is  not  distorted  much  in  smooth  regions,  but 
the  discontinuities  are  spread  over  distance  scales  resolvable  on  a  practical  computational 
mesh.  _  .  ' 

This  approach  leads  to  some  difficulties.  The  width  of  the  transition  layer  representing 
a  discontinuity  in  case  of  a  shock  is  only  few  meshsizes.  However,  it  grows  indefinitely 
in  space  with  the  distance  from  the  initial  boundary  in  case  of  a  contact  discontinuity. 
Another  defect  is  that  when  the  difference  scheme  is  more  than  first  order  accurate, 
this  transition  is  not  monotonic  anymore.  Large  numerical  overshoots  and  undershoots 
occur  from  both  sides  of  the  discontinuity.  Whether  one  needs  to  remove  these  oscillations 
except  for  aesthetic  reasons  seems  to  be  problem  dependent.  For  instance,  one  unpleasant 
consequence  of  these  oscillations  may  be  that  some  quantities  may  attain  values  out  of 
their  physical  domain  (like  negative  density  or  pressure  etc.).  Another  consequence  may 
be  that  in  case  of  a  strong  shock  (when  the  characteristic  field  is  essentially  convergent) 
these  oscillations  may  cause  convergence  to  a  non-physical  solution  or  even  the  loss  of 
stability. 

There  exists  a  whole  family  of  one-dimensional  methods  which  produce  a  higher  order 
solution  in  the  smooth  regions  and  give  a  sharp  resolution  of  discontinuities.  They  are  so- 
called  high  resolution  schemes,  which  are  shown  to  be  total  variation  diminishing  -  TVD 
(see  (11,  12,  16,  21]),  and  lately  essentially  non-oscillatory  schemes  -  ENO  (see  [13]). 
These  methods  are  successfully  used  to  solve  ID  time-dependent  problems.  However, 
even  TVD  schemes  are  shown  to  be  convergent  for  finite  time  only. 

The  situation  with  two-dimensional  methods  is  much  less  favorable.  The  explanation 
is  probably  that  the  physics  of  one-dimensional  flows  is  especially  simple  and  well  under¬ 
stood,  and  easy  to  imitate  by  a  numerical  process.  Also,  it  has  been  shown  in  [9],  that  a 
higher  order  accurate  2D  scheme  cannot  be  TVD.  We  can  conclude  that  the  theoretical 
results  for  the  TVD  schemes  may  not  be  relevant  for  the  higher  order  methods  applied 
to  2D  steady  state  problems,  which  is  our  interest  here. 

One  approach  toward  the  numerical  simulation  of  multidimensional  flows  is  called 
“dimensional  splitting”.  The  multidimensional  differential  operator  is  expressed  as  a 
sequence  of  one-dimensional  operators,  which  are  then  approximated  by  one-dimensional 
difference  operators  (a  review  can  be  found  in  [23]).  This  approach  completely  ignores 
the  multidimensional  nature  of  the  flow  (vorticity,  the  possible  use  of  diagonal  neighbors 
in  the  difference  scheme,  the  infinite  number  of  wave  propagation  directions  instead  of 
only  two  in  ID).  Therefore  it  is  important  to  develop  fully  multidimensional  methods. 

The  key  here  is  probably  to  imitate  the  anisotropic  nature  of  the  multidimensional 
fluid  flow.  This  is  possible  to  achieve  by  using  moving  grids  with  gridlines  aligning 
with  the  stream  direction.  However,  this  simple  idea  turns  out  to  be  very  difficult  to 
implement  even  in  case  of  a  simple  problem.  Our  approach  here  is  to  construct  genuinely 


2D  difference  schemes  which  use  a  fixed  grid,  but  provide  a  separation  between  the 
treatment  of  stream  wise  and  cross-stream  directions.  This  allows  us  to  combine  in  one 
difference  scheme  properties  of  stability  with  good  resolution  of  discontinuities  and  higher 
order  accuracy. 

-  We  show  in  this  work  that,  once  a  discontinuity  is  recognised  in  the  shock  capturing 
solution,  its  location  can  be  recovered  up  to  the  same  order  of  accuracy  as  achieved  by 
the  numerical  solution  in  smooth-  regions.  This  means  that  using  the  shock  capturing 
approach  one  can  obtain  as  much  information  about  a  solution  as  using  the  shock  fitting 
approach/  The  multigrid  solver  apears  to  be  as  efficient  in  converging  the  discontinuity 
location  as  in  converging  the  solution  in  smooth  regions. 

The  next  chapter  presents  a  discussion  about  the  boundary- value  problem  for  the 
2D  non-linear  steady-state  conservation  law,  basic  properties  of  its  solution  and  some 
considerations  regarding  its  discretization.  We  show  that  if  conservative  second  order 
accurate  stable  scheme  is  used,  it  is  possible  to  recover  the  discontinuity  location  from 
the  obtained  numerical  solution  with  second  order  accuracy. 

Chap.3  contains  a  description  a  traditional  multigrid  algorithm  for  elliptic  problems  as ' 
well  as  the  discussion  about  adaptation  of  coarsening  procedure  for  the  non-elliptic  case 
with  disconinuous  solutions,  assuming  that  a  stable  second  order  accurate  discretization 
is  employed  by  the  algorithm. 

Chap.4  is  devoted  to  the  construction  of  such  a  monotonic  discretization  for  the  case 
of  a  linear  constant  coefficient  equation.  The  generalization  of  this  discretization  for  the 
general  non-linear  case  is  presented  in  Chap.5. 

Chap. 6  reports  about  the  numerical  experiments,  where  we  examine  the  accuracy  of 
the  obtained  solutions  (both  in  terms  of  the  solution  error  in  smooth  regions  and  the 
error  in  discontinuity  location  as  well  as  the  perforemance  of  the  algorithm. 

Chap.7  contains  the  summary,  the  efficiency  comparison  with  the  existing  methods 
together  with  the  discussion  of  future  possible  developments. 


Chapter  2 

Conservation  law  and  its  discretization 

A  simple  differential  equation,  but  typical  to  more  complicated  systems  in  fluid  dynamics, 
is  the  2D  nonlinear  steady-state  conservation  law 

-  e&u  +  (/(«)),  +  (g(u))v  =  a(x,  y ),  (2.1) 

where  e  >  0  is  small,  A  denotes  the  laplacian,  /,  g  and  s  are  given  functions. 


2.1  Boundary- value  problem:  solution  properties 

We  shall  discuss  properties  of  the  solution  of  the  boundary-value  problem  for  (2.1)  in 
order  to  get  insight  to  its  discretization. 

2.1.1  Linear  convection  diffusion  equation 

Consider  first  a  linear  version  of  (2.1)  -  a  convection-diffusion  equation 

—  eAu  +  aux  +  buy  =  s(x,  y),  (2.2) 

The  line  whose  tangent  at  every  point  is  determined  by  the  vector  (a(x,y),6(x,y))  is 
called  a  characteristic  line.  The  entire  domain  can  be  covered  by  a  family  of  characteristic 
lines  (characteristics).  We  want  to  single  out  the  time-like  direction  along  these  lines. 
In  order  to  be  consistent  with  the  sign  of  the  elliptic  term  coefficient  the  direction  of 
the  vector  (a,  b)  must  be  the  choice.  We  shall  call  it  the  characteristic  direction  or  the 
stream  direction. 

Consider  the  degenerate  case  ( e  =  0)  of  Eq.(2.2).  This  equation  will  be  hyperbolic 
with  respect  to  the  characteristic  direction  (a,  6).  The  boundary  value  problem  for  this 
equation  is  not  well  posed.  On  the  other  hand,  suppose  we  select  the  part  of  the  boundary, 
at  every  point  of  which  the  vector  (a(x,y),6(x,y))  is  directed  into  the  domain.  Letting 
the  data  on  this  part  of  the  boundary  serve  as  initial  data  for  the  hyperbolic  equation, 
the  obtained  initial-value  problem  for  the  hyperbolic  equation  is  well  posed,  and  its 
solution  has  the  property  that  the  information  from  the  initial  data  propagates  along 
characteristics.  This  means  that  the  value  of  the  solution  at  every  particular  point  of 
each  characteristic  line  depends  only  on  the  initial  value  at  that  point  of  the  boundary 
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where  this  characteristic  line  comes  from  and  on  the  values  of  the  right-hand-side  s  along 
this  characteristic.  Hence,  if  the  initial  data  are  oscillatory,  these  oscillations  will  be 
convected  into  the  domain  along  characteristics  without  any  damping.  We  shall  call  the 
solution  component  which  is  smooth  in  the  streamwise  direction,  but  oscillates  otherwise, 
the  characteristic  component.  If  the  initial  data  are  discontinuous  at  some  point,  the 
solution  of  the  equation  wiU  have  a  discontinuity  (called  a  contact  discontinuity)  along 
the  characteristic  line  which  emanates  from  this  point  of  the  boundary. 

Return  now  to  the  original  elliptic  (viscous)  problem  (2.2).  There  exists  a  layer  of 
width  Q(rrw)  (called  a  boundary  layer),  along  that  part  of  the  boundary  not  taken  as 
the  initial  data  for  the  hyperbolic  problem,  where  the  solution  of  the  elliptic  problem 
may  have  rapid  changes.  The  solution  of  the  original  elliptic  problem  differs  from  the 
solution  of  the  above  described  hyperbolic  problem  only  by  0(e)  in  the  whole  domain, 
except  for  the  boundary  layer.  If  the  solution  of  the  hyperbolic  problem  contains  a 
contact  discontinuity,  the  solution  of  the  original  elliptic  problem  will  have  a  fast  change 
smeared  over  a  layer  of  width  0(r%  ([„[+[E|)^  )>  where  r  is  the  distance  along  characteristic 
from  the  initial  data. 

The  characteristic  components  of  the  solution  will  be  dampened  as  they  are  convected 
along  characteristics.  We  want  to  obtain  a  quantitative  description  of  this  damping.  We 
shall  recall  the  result  for  the  discrete  approximation  of  (2.2)  presented  in  Sec.2  of  [4].  The 
characteristic  component  with  cross-stream  wavelength  rj  will  lose  a  substantial  fraction 
of  its  amplitude  when  reaching  a  certain  distance  r^(rj)  from  the  boundary  into  the 
domain.  This  r*, p{q)  is  called  a  survival  depth  of  the  rj  component  and  as  it  is  shown  in 

S^2  °l  W  rkM  =  0(alV»'h-’),  (2-3) 

where  h  is  the  meshsize,  p  is  the  order  of  approximation,  and  ah  may  reflect  the  “width” 
of  the  stencil,  i.e.  its  diameter  when  projected  on  a  plane  perpendicular  to  the  charac¬ 
teristic  direction. 


2.1.2  Nonlinear  equation 


Consider  now  the  degenerate  case  of  Eq.(2.1).  We  shall  obtain  the  hyperbolic  equation 


(/(«*))*  +  =  (2*4) 

Rewrite  the  equation  (2.4)  in  the  quasilinear  form 

a(u)u9  +  b(u)uv  =  s(x,y),  (2.5) 


where 


a(u)  =  df/du 
b(u)  =  dg/du. 


(2.6) 


As  in  the  linear  case,  characteristics  may  be  introduced  and  the  characteristic  direction 
can  be  determined.  The  part  of  the  boundary  from  which  the  vector  (a,  6)  is  directed 
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into  the  domain  can  be  selected  and  the  boundary  conditions  there  can  be  considered  as 
the  initial  conditions  for  Eq.(2.4). 

Integrate  equation  (2.4)  over  domain  Cl  and  apply  Gauss  theorem 

-  JJQsdxdy ~  JJ^f(u))*^(9^))v\dxdy  =  J^f,g)nd(dCl).^  (2.7) 

Here  dfl  denotes  the  boundary  of  the  domain  fl,  and  ( ,  )n  the  outward  normal  (to  dCl) 
of  the  vector  (,-').  Relation  (2.7)  expresses  conservation:  the  integral  of  the  flux  through 
the  boundary  is  equal  to  the  integral  of  the  sources  inside  the  domain,  u  is  called  a 
generalized  solution  of  Eq.  (2,4)  if  it  satisfies  its  integral  form  (2.7),  i.e.  if  (2.7)  holds 
for  every  smoothly  bounded  domain  Cl. 

The  solution  of  the  initial  value  problem  for  (2.4)  as  well  as  for  (2.7)  may  be  non¬ 
unique  and  discontinuous  even  in  case  of  continuous  initial  data  (see  [10,  22]).  However, 
we  are  interested  only  in  the  solution  which  can  be  obtained  from  the  solution  of  the 
boundary  value  problem  for  Eq.(2.1)  at  the  limit  e  — *■  0.  Such  a  limit  solution  may  be 
discontinuous,  but  is  unique  (see  (22]). 

Consider  again  Eq.(2.5).  The  main  difference  from  the  linear  case  is  that  the  position 
of  characteristics  depends  also  on  the  solution.  Therefore,  characteristics  may  intersect, 
i.e.  information  from  different  points  of  the  initial  data  will  be  brought  along  charac¬ 
teristics  to  the  same  point,  their  intersection  point.  In  such  a  case  discontinuity  of  the 
solution,  which  is  called  a  shock  wave,  will  be  produced.  Only  such  discontinuities  are 
physically  relevant.  In  other  words,  through  every  point  of  the  path  of  a  discontinuity 
in  the  (i,  y)  plane  one  can  draw  two  characteristics,  one  on  each  side  of  the  shock.  The 
discontinuity  will  be  physically  relevant  if  these  two  characteristics  can  be  traced  in  the 
upstream  direction  back  to  the  initial  data.  (Contact  discontinuity  is  a  limit  case  of  such 
a  situation.) 

Introduce  the  following  symmetric  vector- valued  function  of  two  variables 

S(u,t>)  =  [(/(«),  $(u))  -  (/(v), g(v)))  sign(u  -  v).  (2.8) 

Let  T(u)  €  Cl  be  a  set  of  points  where  the  function  u(x,y)  is  discontinuous.  Then  this 
discontinuity  will  be  an  admissible  one  (physically  relevant)  if  the  following  inequality 
(the  entropy  condition)  is  satisfied 

(S(u+,  c),  v )  <  (S(tt~,  c),  y),  (2.9) 

where  c  is  an  arbitrary  constant,  v  is  the  normal  to  T(u),  u+  is  the  limit  value  of  the 
solution  from  the  side  of  discontinuity  where  v  is  directed  to,  u~  the  limit  value  of  the 
solution  from  the  opposite  side  of  the  discontinuity  ([22]). 

The  generalized  solution  of  (2.4),  which  satisfies  (2.9)  is  unique  and  it  coincides  with 
the  limit  solution  (when  e  — ►  0)  of  the  boundary  value  problem  for  (2.1)  ([22]), 
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2.2  Discretization  of  conservation  law 

Our  approach  to  the  discretization  of  (2.4)  will  be  guided  by  its  integral  form  (2.7), 
because  it  holds  also  accross  discontinuities. 

Consider  a  computational  grid  with,  a  meshsize  A  in  both  x  and  y  directions v  which 
covers  domain  Cl.  Integrate  Eq.(2.4)  over  computational  cell  A,  surrounding  gridpoint 

(*,i)  (se®  Fig*  2.1)  ' 

L  (MnWCuU  //  ,<tzdy.  (2.10) 

The  integral  in  the  right-hand-side  of  (2.10)  can  be  splitted  into  four  parts  -  one  along 
each  segment  of  the  computational  cell  boundary.  Since  we  are  interested  in  second  order 
accurate  solutions,  it  will  be  sufficiently  accurate  to  approximate  these  integrals  using 
the  mid-point  rule.  The  2D  mid-point  rule  can  be  used  to  approximate  the  left-hand-side 
of  (2.10)  as  well. 

Then  the  discrete  form  of  (2.10)  (we  shall  call  it  the  balance  equation)  can  be  written 
as  follows: 

^  +  9*J+$  ~  3ij-\)  —  h*  3iJ  (2-11) 

The  difference  scheme  constructed  this  way  is  called  conservative  in  that  it  has  a  property 
analogous  to  (2.7):  in  the  sum  of  discrete  equations  over  all  the  grid  points  only  numerical 
fluxes  through  the  boundary  remain  and  all  other  fluxes  cancel  each  other.  The  method 
of  calculating  the  numerical  fluxes  in  (2.11)  will  be  presented  in  Chapters  4,5. 


2.3  Recovering  of  discontinuity  location 

As  well  known,  the  basic  advantage  of  conservative  schemes  is  that  when  the  meshsize 
A  — ►  0  the  discontinuities  will  converge  to  their  correct  location.  However,  the  question 
of  the  order  of  convergence  has  remained  open.  Moreover,  it  was  not  clear  what  should 
be  understood  by  a  discontinuity  location  in  case  of  a  discrete  solution.  We  define  here 
a  discontinuity  location  for  a  shock  capturing  solution  and  we  show  that  A2  convergence 
in  it  can  be  achieved,  if  there  is  A 2  convergence  in  the  smooth  regions. 

2.3*1  Integral  relation 

Suppose  the  exact  physically  relevant  solution  of  (2.4)  contains  a  discontinuity  with  a 
path  d  (see  Fig.  2.2). 

Suppose  we  also  have  a  numerical  solution  of  (2.4)  obtained  by  means  of  a  certain 
conservative  finite  difference  scheme  (2.11).  The  discontinuity  will  be  represented  in  the 
numerical  solution  by  a  transition  layer  (monotonous  or  containing  oscillations).  Assume 
that  the  numerical  fluxes  approximate  the  exact  fluxes  at  the  corresponding  points  with 
accuracy  0(h2)  in  the  smooth  regions,  away  from  the  influence  of  the  transition  layer. 


•n 
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Figure  2.2:  Domain  and  discontinuity  path. 

Draw  a  segment  LR  across  the  discontinuity.  Suppose  that  this  segment  goes  along 
computational  cell  boundaries  and  that  there  is  already  a  pointwise  h 2  convergence  to 
the  exact  solution  at  the  point  L  and  R.  We  want  to  connect  each  of  these  two  points 
with  the  boundary  by  a  path  which  will  go  along  computational  cell  boundary  and  will 
belong  to  the  region  with  h 7  convergence.  Keeping  these  lines  close  to  characteristics 
going  from  L  and  R  back  to  the  boundary  (Fig.  2.2)  seems  to  be  a  good  choice,  because 
of  the  physical  relevance  of  the  solution.  The  integration  of  (2.4)  over  the  subdomain 
bounded  by  BlLRBr  using  Gauss  theorem  will  give 

^(/,y)n  d{dw)  =  Jjusdx  dy-  (2.12) 

Split  the  boundary  integral  into  four  parts: 

/  if,g)nd(du>)  = 

JOw 

L  a  (f>9)nd(du)  +  /  (f,g)nd{du)  + 

JBrBl  JBlL 

(  {f,g)»d{d u/)  +  f  B(f,g)nd(du>).  (2.13) 

JRBr  Jlr 

By  (2.11)  the  summation  of  the  balance  equations  over  the  computational  cells  in¬ 
cluded  in  u  is 

h  ”  fi-y  +  9iJ+i  -  g*j-\)  =  h'2'52  Ki-  (2.14) 

Notice  that  in  the  last  summation  all  numerical  fluxes  cancel  each  other  except  for  those 
through  the  boundary  7  of  subdomain  u.  Therefore,  we  denote 

9u  w 


(2.15) 
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Figure  2.3:  Discontinuity  location  recovering. 


Since  we  can  asstime 

JJ  s  dxdy  —  hi2  ^2  Sij  +  0(h2),  (2.16) 

we  can  deduce  from  (2.15),(2.14)  and  (2.12)  that 

l  (f,9)nd(du)  =  hJ2(f\gh)n  +  0(h2).  (2.17) 

J9v  aw 

But  since  there  is  0(h2)  pointwise  convergence  along  the  lines  BrBj^BiL  and  RBr,  on 
the  remaining  line  LR  we  must  have 

I  Af,  g)n  d(d u)  -  h  £(/\  gh)n  =  0{h 2),  (2.18) 

JLR  LR 

despite  the  lack  of  pointwise  second  order  convergence  along  this  segment  LR. 

2.3.2  Implementation 

We  can  assume  without  loss  of  generality  that  the  segment  LR  is  perpendicular  to  the 
y-direction.  Then,  the  flux  normal  to  this  segment  will  be  just  g( u). 

The  numerical  fluxes  values  gi  and  gR  are  assumed  to  approximate  the  fluxes  of  the 
differential  solution  at  the  respective  points  with  an  accuracy  of  h 2.  The  values  at  points 
between  L  and  R  do  not  in  general  approximate  the  differential  solution.  Their  only  role 
is  to  indicate  the  discontinuity  location. 

Let  us  reconstruct  fluxes  of  a  discontinuous  solution  from  the  numerical  solution  we 
have.  For  this  purpose  we  have  to  substitute  the  values  in  the  transition  layer  by  values, 
which  approximate  the  differential  solution.  Therefore,  we  have  to  extrapolate  them 
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from  the  grid  points  to  ne  left  of  L  and  to  the  right  of  R  into  the  segment  LR  (see  Fig. 
2.3).  The  values,  extrapolated  from  the  left  will  approximate  fluxes  of  the  differential 
solution  on  the  left  side  of  the  discontinuity,  and  values  extrapolated  from  the  right  will 
aproximate  fluxes  from  the  right  side  of  the  discontinuity.  We  just  have  to  find  the  proper 
location  of  the  discontinuity  (on  this  segment  LR).  Denote  the  reconstructed  fluxes  g^. 
We  shall  locate  the  discontinuity  in  such  a  way,  that  the  following  relation  will  hold 

7  (2.19) 

-  ~  JLR  -  LR  -  - 

This  relation  defines  the  discontinuity  location  to  within  0(h2)  accuracy.  Indeed,  suppose 
first  that  the  width  of  the  transition  layer  behaves  like  0(h).  This  means  that  the  number 
of  grid  points  involved  in  the  transition  layer  remains  the  same  on  each  grid.  Then,  in 
order  to  obtain  A 2  accuracy  in  the  discontinuity  location  it  is  enough  to  use  just  an 
extrapolation  by  constants  from  points  L  and  R.  In  case  the  width  of  the  transition 
layer  behave  like  O(A^),  it  is  enough  to  use  a  linear  extrapolation. 

The  argument  presented  here  has  important  implications  for  multigrid  methods. 
These  implications  are  explained  in  Chap.3. 


Chapter  3 

Multigrid  methods 


In  this  chapter  we  first  give  a  brief  description  of  the  multigrid  solver  for  an  elliptic 
problem.  Then  we  show  how  to  adapt  a  multigrid  algorithm  for  a  non-elliptic  problem. 

3.1  Elliptic  case 

3.1.1  Relaxation 

Consider  a  difference  scheme 

LV  =  sk  (3.1) 

with  certain  boundary  conditions,  approximating  a  boundary-value  problem  for  an  ellip¬ 
tic  differential  equation  Lu  =  s. 

Suppose  these  difference  equations  are  being  solved  by  a  certain  relaxation  performed 
in  a  certain  ordering.  The  error  after  n  relaxation  sweeps  is 

=  (3.2) 

where  u*  is  the  current  solution  approximation. 

It  can  be  observed  in  numerical  experiments  that  convergence  of  the  relaxation  is  fast 
in  the  beginning,  but  becomes  very  slow  after  few  sweeps.  This  is  because  the  relaxation 
appears  to  be  very  efficient  in  reducing  of  the  non-smooth  error  components.  When  the 
error  is  smooth  the  convergence  is  slow.  That  is  relaxation  smoothes  the  error.  However, 
the  smooth  error  can  be  approximated  on  a  coarser  grid.  This  is  the  main  idea  of  the 
multigrid  methods. 

3.1.2  Coarse  grid  correction 

Assume  the  Eq.(3.1)  to  be  linear.  Let  uk  be  the  approximation  obtained  by  a  few 
relaxation  sweeps  on  the  fine  grid.  The  residuals  of  the  fine  grid  equations  are  then  given 
by 

rk  s  ak  -  Lkuk  (3.3) 


The  error  v*  =  u*  —  uk  will  satisfy 


LV  =  r*. 


(3.4) 
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Since  vk  is  smooth  it  can  be  approximated  by  a  coarse  grid  function  vH  satisfying  the 
equation 

Lhvh  =  th  ,  (3.5) 

where  LH  is  a  grid  H  difference  approximation  to  the  differential  operator  and 

-  .  -  '  rH  =  /f  r\  .  *  (3.6) 

where  iff  is  a  fine-to-coarse  transfer  operator  (injection  or  a  certain  type,  of  weighting). 

After  obtaining  a  certain  approximate  solution  vH  of  (3.5)  we  can  use  it  to  correct 
the  fine  grid  solution  in  the  following  way 

uA  <-  uk  +  J&t;\  (3.7) 

where  Ig  is  an  interpolation  operator. 

The  process  of  calculating  rk,  transfering  it  to  the  coarse  grid,  solving  the  coarse  grid 
equations  for  v11  and  interpolating  the  result  and  adding  it  to  the  fine  grid  approximation 
is  called  “coarse  grid  correction”. 


3.1.3  Full  approximation  scheme 

We  have  described  a  way  to  construct  coarse  grid  equations  for  the  error  vk  in  the 
fine  grid  approximation  u*  to  the  solution  of  (3.1).  It  is  called  the  Correction  Scheme 
(CS).  Another  way  to  do  it  is  called  the  Full  Aproximation  Scheme  (FAS).  Coarse  grid 
equations  are  written  in  terms  of  a  different  function:  instead  of  vH  we  use 

uH  =  uk  +  vH.  (3.8) 

This  coarse  grid  function  approximates  the  full  solution  on  the  coarse  grid.  The  equation 
for  it  is 

Lhmh  =  s*  (3.9) 

where 

SH  =  LHI?uk  +  If!  r\  (3.10) 

After  solving  this  equation  approximately,  we  use  vF  to  correct  the  fine  grid  approxima¬ 
tion  in  the  following  way 

uk  <-  uk  +  /&( uH  -  i*uk).  (3.11) 

FAS  is  used  in  case  of  nonlinear  problem  or  when  local  refinement  is  needed. 

Both  CS  and  FAS  schemes  can  be  applied  recursively.  Such  a  solution  method  will 
be  very  efficient  not  only  because  coarser  grid  consist  of  less  grid  points.  The  error 
components,  which  are  smooth  on  the  finer  grid,  will  look  “less”  smooth  on  coarser  grids 
and  therefore  can  be  efficiently  reduced  there  by  relaxation. 
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3.1.4  Multigrid  cycle 

Assume  we  have  a  sequence  of  grids  hi  >  hi  >  ...  >  where  h\{  =  h  and  hk  =  2hk+i. 
The  k  grid  equation  is  written  as 

Lkuk  =  sk.  (3.12) 

For  k  =  M:  Lk  =  Lk  and  sk  =  ah.  In  the  Case  of  the  Correction  Scheme 

**■*  =  ’■  ■  (3.13) 

and  the  initial  approximation  for  it*"1  is  zero.  For  the  Full  Aproximation  Scheme 

3*-1  =  Lk-xIkk~'uk  +  Ikk'\sk  -  Lkuk )  (3.14) 

and  the  initial  approximation  for  u*-1  is 

Given  an  approximate  solution  uk  of  the  Eq.(3.12).  The  multigrid  cycle 

u*  —  MG(k,uk,sk)  (3.15) 

is  defined  recursively  as  follows: 

•  If  k  =  1,  solve  (3.12)  by  relaxations. 

•  Otherwise: 

—  Perform  v\  relaxation  sweeps  on  (3.12)  resulting  in  a  new  approximation. 

-  Transfer  the  problem  to  the  coarser  grid  1  —  1  and  perform  7  successive  cycles 

uk~l  «-  M G[k  - 1,  uk-\sk-1).  (3.16) 

—  Interpolate  the  correction  to  grid  k  and  add  it  to  the  approximate  solution 
there. 

-  Perform  vj  relaxation  sweeps  on  grid  k  resulting  in  u*  of  (3.15). 

When  7  =  1  the  corresponding  multigrid  cycle  is  called  V  cycle  or  V(vi,  v-i).  If  7  =  2 
-  it  is  called  W  cycle  or  W(  17,  i/j). 


3.1.5  Full  multigrid  algorithm 

The  first  approximation  for  the  multigrid  cycle  on  a  certain  level  can  be  obtained  by 
interpolating  a  solution  from  the  previous  coarser  level,  which  itself  has  been  calculated 
by  multigrid  cycles  MG  on  that  level  and  so  on  recursively.  Such  an  algorithm  is  called 
a  Full  Multigrid  Algorithm  ( FMG ).  We  shall  denote  by  FM G (N,  Nm  ,C,M)  a  certain 
FMG  algorithm, where  M  is  the  finest  level,  Nm  -  the  number  of  multigrid  cycles  per¬ 
formed  in  this  finest  level,  N  -  number  of  cycles  performed  on  the  intermediate  levels 
and  C  is  the  type  of  a  multigrid  cycle  employed  (V(uu  47)  or  W{v 1,1^3)).  In  case  we 
want  only  to  clarify  how  many  MG  cycles  are  done  on  the  finest  level  we  shall  use  the 
notation  Nm  —  FMG. 
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3.2  Non-elliptic  case 

The  first  task  in  this  case  is  to  construct  a  stable  discretization,  such  that  a  certain  non- 
expensive  local  relaxation  can  be  applied  with  it  and  will  have  good  smoothing  properties. 
It  is  sufficient  for  this  purpose  if  the  discrete  equations  are  ^-elliptic  (see  [2,  3]),  which 
can  be -achieved  by  adding  artificial  viscosity  with  the  coefficient  proportional  to  h. 
Of  course,  this  restricts  us  to  the  first  order  accuracy,  hut  once  it  is -done,  the  usual 
multigrid  algorithms  can  be  efficiently  applied  in  case  of  a  smooth  solution.  However, 
the  discontinuous  cases  and  the  possibility  to  obtain  second  order  accuracy  require  further 
studies.  We  require  from  our  discretization  to  be  stable,  second  order  accurate  and  to 
provide  a  good  resolution  of  discontinuities  (representing  them  by  thin  and  monotonic 
transition  layers).  Chapters  4,5  are  devoted  to  the  construction  of  such  discretization 
schemes. 

Since  the  difference  schemes,  which  will  be  developed  in  Chapters  4,5  are  nonlinear 
even  in  case  of  a  linear  differential  equation,  the  Fhll  Approximation  Scheme  should 
be  used.  We  have  shown  in  Chap.2  what  should  be  understood  by  a  discontinuity 
location  when  a  shock  capturing  approach  is  used.  The  question  is  now  how  to  perform 
coarsening  in  the  multigrid  process  in  order  to  achieve  the  same  efficiency  in  converging 
the  discontinuity  to  its  correct  location  as  in  converging  the  solution  in  the  smooth 
regions. 

We  have  demonstrated  that  the  conservation  property  of  the  difference  scheme  is  of 
crucial  importance  for  obtaining  the  correct  discontinuity  location.  Therefore,  in  order 
to  obtain  a  proper  coarse  grid  correction  for  the  discontinuity  location  a  conservative 
residual  transfer  should  also  be  used.  This  can  be  just  the  usual  Full- Weighting. 

Another  question  is  what  correction  interpolation  should  be  used  for  this  purpose. 
The  important  consideration  here  is  preserving  the  flux  correction  integral  along  coarse 
grid  computational  cells  boundaries.  Since  this  correction  is  small  (when  an  FMG 
algorithm  is  used  it  is  0(h2)  in  smooth  regions  and  0(h)  in  the  neighborhood  of  a 
discontinuity),  it  is  enough  to  preserve  integral  of  the  solution  correction  itself.  This  is 
perfectly  done  by  the  usual  bilinear  interpolation,  which  was  shown  already  to  perform 
well  in  smooth  regions.  Therefore,  the  conclusion  is  that  for  obtaining  a  full  multigrid 
efficiency  for  converging  the  discontinuity  location  the  usual  coarsening  techniques  can 
be  used.  In  other  words,  the  coarse  grid  correction  is  perfectly  capable  of  moving  the 
discontinuity,  commensurably  with  the  solution  changes  it  affects.  There  is  no  need  to 
freeze  the  discontinuity  location  before  going  to  coarse  grids  and  performing  a  special 
procedure  to  update  it  after  coming  back  to  the  fine  grid  (like  it  was  suggested  in  [5]). 


Chapter  4 


Discretization:  Linear  case 


The  problem  we  shall  deal  with  now  is  how  to  construct  numerical  fluxes 

The  general  nonlinear  case  will  be  considered  in  the  next  chapter.  In  the  present 
chapter  we  shall  study  some  existing  approaches  as  well  as  our  own  on  the  case  of  the 
simple  linear  constant  coefficient  equation 

—  eAu  +  au*  +  buy  =  s,  (4.1) 

where  s  =  s(x,  y).  Without  loss  of  generality  we  can  assume 

&  >  a  >  0,  6  >  0.  (4.2) 

Since  the  fluxes  fi+^  and  /,_^j  are  constructed  in  the  same  way  (as  well  as 
and  we  3hall  give  formulas  for  /,_^j  and  9ij-$  only- 


4.1  Some  existing  methods 


4.1.1  Central  differencing 

The  most  straightforward  approach  is  the  following  “central’’  fluxes: 

=  i&Ki  +  “iJ-i)-  '  ’ 

When  (4.3)  is  substituted  into  the  balance  equation  (2.11),  it  will  give  the  standard 
central  5-point  “star”  second  order  accurate  approximation  to  the  equation  (4.1) 

I  -  «-*■<  +  t«W  -  -vr! )  ,  (4.4) 

This  scheme  does  not  have  a  good  measure  of  ellipticity  (see  [3]).  There  exist  certain  high- 
frequency  components  which  can  be  present  in  the  error,  but  do  not  express  themselves 
in  residuals. 


4.1  Some  existing  method* 
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4.1.2  Upstream  differencing 


Let  us  add  additional  terms  to  the  fluxes  defined  by  (4.3)  to  produce  “upstream”  fluxes: 


/V*J  =  ~  «*-v)  =  a*-U  ,,  ,, 

9uij-±  =  9cij-±  .  -  *  ’ 

These  additional  terms  correspond  to  the  following  anisotropic  artificial  viscosity 

-  \K(a**)*  +  (K)v)>  -  ;  '  (4.6) 


which  has  the  same  sign  as  the  physical  viscosity.  When  fluxes  (4.5)  are  substituted 
into  the  balance  equation,  we  obtain  a  first  order  accurate  upstream  scheme,  based  on  a 
3-point  stencil  which  can  also  be  written  in  the  form 


=  7+6u‘"u  +  TTb^-" 


(4.7) 


Since  the  coefficients  in  the  right-hand  side  of  this  equation  are  positive  (i.e.  this  differ¬ 
ence  operator  is  of  the  positive  type),  the  following  relation  will  hold: 


min  («,-_!  j,  liij-i)  <  <  max  (u,_u,  Uij.%).  (4.8) 

In  other  words,  a  certain  maximum  principle  holds  for  the  solution  of  (4.5).  We  can 
summarize  that  the  addition  of  the  artificial  viscosity  to  the  central  second  order  scheme 
creates  a  stable  upstream  first  order  scheme  which  produces  a  monotonic  solution.  The 
characteristic  components  (i.e.,  components,  which  are  smooth  in  the  characteristic  di¬ 
rection  and  oscillating  otherwise)  and  discontinuities  propagating  in  an  oblique  direction 
will  be  smeared  significantly  in  the  solution,  but  those  propagating  in  the  grid  direction 
(i.e.,  when  a  =  0)  will  be  resolved  perfectly. 


4.1.3  Upstream  second  order  scheme 


We  would  like  to  construct  a  second  order  accurate  scheme,  which  will  maintain  the 
stability  property  of  (4.5).  We  can  approximate  derivatives  in  each  direction  using  second 
order  accurate  3-point  one-sided  approximation.  The  difference  equation  will  be 


3tiij  -  4ui_u  +  ,  t3uij  —  4u,j_l  -f-  u,j_2 

1  2  h  +  2  h 


’•j- 


(4.9) 


In  terms  of  numerical  fluxes  it  can  be  written 

=  /  Vf  j  + 

~  +  I&Km  -  vij-*), 


(4.10) 


The  role  of  the  additional  terms  is  to  compensate  for  the  loss  of  accuracy  due  to  the  arti¬ 
ficial  viscosity  (diffusion).  Therefore,  they  are  often  called  antidiffusive  fluxes.  However, 
when  applied  in  the  neighborhood  of  a  discontinuity,  this  scheme  will  produce  spurious 
oscillations  in  the  solution. 


Chapter  4  —  Discretization;  Linear  case 
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4.1.4  Remark  on  overshoots 

Consider  a  first  order  scheme.  Its  truncation  error  is  dominated  by  its  artificial  viscosity, 
which  behaves  like  a  physical  one.  Therefore,  the  truncation  error  causes  at  least  as 
much  second  order  dissipation  at  every  point  as  the  .physical  viscosity.  As  a  result, 
the  difference  solution  has  a  mono  tonicity  property  analogous  to  one  of  the  differential 
solution  and  overshoots  do  not  appear.  In  case  of  the  higher  order  difference  scheme, 
its  truncation  error  is  dominated  by  higher  order  derivatives  and  does  not  behave  like  a 
physical  viscosity  anymore.  It  may  even  attain  values  which  are  opposite  in  sign  to  that 
of  the  second  order  dissipative  term  at  some  grid  points.  This  means  that  some  non¬ 
physical  amounts  of  preserved  quantity  (mass,  momentum  or  energy)  are  injected  at  these 
points.  In  the  neighborhood  of  discontinuities,  where  derivatives  become  la'ge,  these 
local  violations  of  the  conservation  law  cause  appearance  of  large  overshoots.  However, 
if  the  difference  scheme  is  conservative,  the  conservation  law  still  holds  globally. 

The  usual  way  to  overcome  this  difficulty  is  to  multiply  antidiffusive  fluxes  by  a  certain 
quantity,  which  is  called  a  limiter.  The  concrete  values  of  the  limiter  at  each  grid  point 
should  be  chosen  from  considerations  of  the  second  order  accuracy  and  monotonicity.  It 
is  closed  to  1  in  the  smooth  regions,  not  distorting  the  order  of  accuracy  of  the  original 
scheme.  However,  it  becomes  different  from  1  in  the  regions  of  changing  gradients: 
introducing  first  order  artificial  viscosity  needed  to  damp  the  oscillations  if  smaller  than 
1,  or  sharpening  transition  layers  representing  discontinuities  (artificial  compression) 
when  larger  than  I.  This  approach  appeared  to  be  very  successful  for  ID  problems  (the 
high  resolution  schemes  are  based  on  this  principle,  see  [11,  12,  14,  16]).  However,  a 
straightforward  extension  for  2D  (dimensional  splitting)  has  some  flaws.  It  leads  to  wide 
stencils  and  may  require  complicated  block-relaxation  process  in  order  to  maintain  its 
stability. 


4.1.5 

Approach  of  Spekreijse 

Let  us  multiply  the  antidiffusive  fluxes  in  (4.10)  by  a  certain  limiter  ip(R),  producing: 

/\-$j  =  A-ij  +  -  w-2j) 

9*ij-±  =  9Uij- 1  +  iWQ'ij- 1  -  «ij- 2). 

(4.U) 

where 

RT^xj- 

«t-Ij  -  «t-2J 

(4.12) 

Q'ij.1  =  — ■ *— ■ J—. 

J  3  Uij-l  - 

(4.13) 

Choose  ip  to  be  the  Van  Albada’s  limiter 

.  &  +  R 

MR)=  IP +  V 

(4.14) 

4.2  2D  approach:  Homogeneous  equation 
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It  is  shown  in  [20]  that  the  solution  of  (4.11)  will  have  the  same  monotonic  property 
(4.8)  as  the  solution  of  (4.5),  and  will  be  second  order  accurate  in  the  smooth  regions. 

However,  a  pointwise  relaxation  will  be  unstable  when  applied  to  this  scheme.  The 
only  way  to  maintain  the  stability  is  tc  use  more  expensive  block  relaxation  (4-points 
.relaxation,  [20]).  This  scheme  can  be  regarded  as  a  typical  example  of  the  dimensional 
splitting  approach.  .  -  •  '  ’ 


4.2  2D  approach:  Homogeneous  equation 

Consider  the  Eq.(4.1)  with  s(x,y)  =  0. 


4.2.1  2D  scheme 

Our  approach  is  to  construct  upstream  antidiffusive  fluxes  using  certain  difference  ap¬ 
proximations  to  the  original  equation  (4.1)  itself.  It  leads  to  a  “compact”  stencil,  which 
will  use  diagonal  grid  points  instead  of  the  points  which  are  distance  2 h  from  the  central 
one.  The  difference  scheme  based  on  such  a  stencil  will  have  a  smaller  truncation  error 
and  will  be  applicable  near  the  boundaries  as  welL  Another  advantage  is  the  possibility 
to  separate  between  the  treatment  of  the  cross-stream  and  streamwise  directions.  This 
allows  us  to  add  the  high  resolution  in  the  cross-stream  direction  only,  but  to  maintain 
a  good  measure  of  eUiptidty  in  the  streamwise  direction.  As  a  result  we  can  obtain  a 
discrete  scheme  which  is  able  to  resolve  well  characteristic  components  (see  Chap.2)  and 
discontinuities  together  with  good  stability  properties. 

Define  the  following  2D  compact  scheme 

f  i— lj 

=  ah  _  2°(u*v-i  -  «.•- i.i-») 

Lemma  4.1  Scheme  (4.15)  is  second  order  accurate . 

Proof:  We  want  to  show  that 

+iJ  -  =»(<■“.) +0(h3) 

J+i-s’Vi  =  MK) + o(h’). 

Rewrite  (4.15)  as 

*  /*(.  -  m-w)  +  w-0) 

**  4*u-)  “  j(*(“U  ~  «u-l)  +  -  “i-W-l)) 


(4.16) 


(4.17) 
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Chapter  4  -  Discretization:  Linear  case 


(4.19) 


Note,  that 

f2D^j-f2D<-u  = 

|(aKwj  -  Uij)  +  + 

|(aKi  - ;«i-w)  -  +  _b(ui-ij  Ut-ij-i))  (4.18) 

Since  the  centralscheme  (4.3)  is  second  order  accurate,  we  can  conclude  that  . 

f2D.  ,  .  _  f*D.  .  .  _ 

/  «+$j  J  t-\j  — 

h(aux)  -  ~ h2(aut  +  buy)s  +  0(h 3),  (4.19) 

or  taking  the  equation  (4.1)  into  account 

f2D  _  f2D  .  _ 

h(aux)  +  0(kz),  (4.20) 

Similarly 

ff3Dij+i  ~  97Dij-l  =  MK)  -  \h2(«u*  +  K)v  +  0(h3)  = 

MK)  +  0(h3).  (4.21) 

Remark  4.1  There  ia  no  first  order  artificial  viscosity  used  by  the  scheme  (4.15). 


(4.20) 


When  (4.15)  are  substituted  into  the  balance  equation  they  give: 

b  —  a  a  — b  . .  _ 

u<j  =  Ui-ij-i  +  7g«i-w-  (4-22) 

Note  that  the  diagonal  grid  point  value  Uj_u_i  participates  in  the  difference  equation 
instead  of  the  values  u,_2j-  and  u,j_ j  in  (4.11).  This  scheme  is  second  order  accurate, 
but  not  of  the  positive  type,  because  the  coefficient  of  Ui-ij  is  negative. 


Remark  4.2  Another  significant  defect  of  this  scheme  is  that  in  the  case  a  -*■  0  it  leads 
to  the  difference  equation 

=  tti-ij-i  +  (Uij_!  -  ut_ij)  (4.23) 

which  is  based  on  the  wide  stencil  However,  it  is  natural  to  expect  the  following  equation 
in  this  case 

«ij  =  Ui-ij-i-  (4.24) 

This  defect  creates  an  obvious  difficulty  for  extending  this  approach  for  the  case  of  the 
variable  coefficient  equation  and  for  the  non-linear  case  as  well.  This  is  because  it  may 
lead  to  discontinuous  numerical  fluxes. 


4.2.2  First  order  N  scheme 


Modify  the  previous  scheme: 

=  aui-ij  ~  jniin(o,  6)(u,_lo-  -  u,_w_i) 
.=  Kj-1- -  u«-u-i) 

In  our  .case  (4.2)  j  (4.25)  can' be  rewritten  as  . 

Ai-i  “  Ki-i  “  W«.Vi-i  - 

When  substituted  into  the  balance  equation  it  will  give 


a  b  —  a 

-u,_ij_x  H — — 


(4.25) 

5 

ir 

(4.26) 

(4.27) 


This  scheme  is  already  of  positive  type,  however  only  first  order  accurate.  This  type  of 
scheme  was  actually  suggested  in  {6,  7],  however,  the  present  formulation  is  much  simpler 
(by  avoiding  rotation  transformation  of  the  equation)  and  is  obviously  conservative. 
We  shall  call  it  N  scheme  (“narrow”  scheme).  Its  solution  will  satisfy  the  following 
monotonidty  property. 


min  (uij-i,  Ui-u-i)  <  t*y  < 


(4.28) 


Notice,  that  this  scheme  will  perfectly  resolve  discontinuities  that  are  aligned  with  a 
diagonal  direction  as  well  as  along  grid  lines.  It  can  be  shown  that  the  cross-stream 
viscosity  coefficient  of  this  scheme  is  at  least  3.64  times  smaller  (see  [18])  than  that 
of  (4.5).  Still,  the  resolution  of  an  oblique  discontinuity  by  such  a  scheme  has  to  be 
improved. 


4.2.3  S  scheme 

We  would  like  to  correct  the  N  scheme  to  be  second  order  accurate  in  a  way  similar  to 
(4.11),  retaining  the  monotonic  property  (4.28).  This  can  be  done  by  the  following  S 
scheme: 


Ah 


where 


*  -  ui-w-i) 

~  9° 

d  u— i) 

Kut-u  -  Vi-ij-i)  ' 


(4.29) 


(4.30) 


The  question  is  what  conditions  have  to  be  imposed  on  the  limiter  function  in  order 
to  ensure  monotonidty  and  second  order  accuracy  of  the  S  scheme.  The  two  following 
IwnwiM  and  the  remark  answer  this  question  and  show  that  all  the  limiters  used  in  ID 
problems  and  reviewed  in  [21]  are  good  for  this  purpose. 


Lemma  4.2  If  the  limiter  $  =  rj>(R)  satisfies  the  following  inequality 


0  <  ^  <  2,  'f{R)  <  2, 

then  the  monotonicity  property  (4.28)  holds  for  the  S  scheme. 


(4.31) 


Proof:  Rewrite  (4.30)  as 


1  a. 


U»-  1J  —  ~"n  1  u«-lj-l)' 

Then  the  correction  for  f°_^j  in  (4.29)  can  be  replaced  by 


(4.32) 


1  ,ln  WL  x,  1  a(d  -  a) ,  N 

-  tti-U-i)  =  ~2~R~~f- - ~  (4-33) 

When  the  numerical  fluxes  (defined  by  (4.29))  and  the  correction  for  the  flux  f°_^j  as 
calculated  using  (4.33)  are  substituted  into  the  balance  equation,  we  obtain  the  following 
relation: 


‘“a**.- 


+(6  -  o)(i  -  WRi+ij) + *Tr#i>  (“K.  ~  “m-0  =  »• 


(4.34) 


It  is  easy  to  see  that  condition  (4.31)  ensures  the  positiveness  of  the  coefficients  in  (4.34), 
and  therefore,  monotonidty  of  the  S  scheme. 

Remark  4.3  It  is  easy  to  see  from  (4.34)  that  the  S  scheme  in  case  a  — ►  0  leads  to  the 
difference  equation  based  upon  the  narrow  stencil 


=  “o-i- 


(4.35) 


This  means  that  the  S  scheme  does  not  suffer  from  the  same  defect  as  the  2D  scheme. 


Lemma  4.3  lf^>  —  €  C2  and 


t£(1)  =  1, 


(4.36) 


then  the  S  scheme  is  second  order  accurate. 


Proof:  In  order  to  prove  this  lemma  it  is  enough  to  show  that 

/44j  -  =  *(<■“.) + Of*3). 


(4.37) 


4.2  3D  approach:  Homogeneous  equation 
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f??y  -  fs<-y  =  |(1  -  *«-}„•))(»  -  «■)(<«-«  -  (4-38) 

Taking  (4.36)  into  account: 

1  -WRi-y)  =  -M^-y  - 1)  +  0(k%  .  ;  (4.39)  . 

and  hence  -  -  ' 

\fSy^  /V*  *  "  1)(6  “WO  +0(**) '  •  (iM1  : 

=  MD 1 — (Hui-lj  —  —  Ui-U-l))  +  0(h3). 

Therefore 

(/#„  ~  /V*j)  -  (/£&, „•  -  /%„)  =  0(h3)  (4.41) 

or 

z?+i,  -  z5.-), = /&,  -  z“H3 + <4-42) 

which  proves  (4.37). 

Remark  4.4  Suppose  that  the  limiter  0(72)  is  twice  continuously  differentiable  only  in 
the  neighborhood  of  72  =  1  and  Lipschitz  continuous  otherwise.  It  is  dear  from  the 
proof  of  the  previous  lemma  that  the  S  scheme  employing  such  a  limiter  is  also  second 
order  accurate.  Moreover,  even  if  0*  is  discontinuous  at  72  as  1,  the  S  scheme  will  be 
still  second  order  accurate  in  L\  norm.  This  is  because  the  order  of  approximation  will 
deteriorate  to  the  first  due  to  discontinuities  in  the  derivative  of  the  limiter  in  sudh  a 
case  only  in  computational  cells  located  along  isolated  characteristic  lines.  Therefore, 
these  computational  cells  will  cover  only  0(h)  part  of  the  domain. 

Remark  4.5  Note  that  a  limiter  satisfying  both  (4.31)  and  (4.36)  cannot  be  a  linear 
function.  Therefore,  a  monotone  higher  order  accurate  scheme  will  be  nonlinear  even  in 
case  of  a  linear  equation. 


4.2.4  Examples  of  limiters 

The  construction  of  a  limiter  is  not  unique.  Therefore,  several  different  limiters  were 
proposed  in  the  ID  framework  (see  [21])  and  can  serve  our  purpose.  The  difference 
between  different  limiters  is  only  in  the  amount  of  artificial  viscosity  they  may  add  to 
the  scheme  in  order  to  damp  oscillations  and  in  the  amount  of  artifidal  compression  they 
may  or  may  not  add  to  the  scheme. 

We  shall  give  here  some  examples  of  limiters. 

Example  4.1  Van  Leer  limiter 

^”1  +  1721  *4,43* 


.  9 
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Chapter  4  -  D  acre  tizat  ion:  Linear  case 


Note  that  xpvL  >  0  is  a  monotone  increasing  function 


.  f  0,  if  R  < 


Example  4.2  A  class  of  limiters  defined  by 

-0*  =  max(0rmin(^/2,  l),min(i2,^)), 

where  1  <  <f>  <  2. 


(4.44) 

(4.45) 


Note  that  Vv  is  a  monotone  increasing  function.  When  <p  =  2,  tp+(R)  corresponds  to  the 
Roe  highly  compressive  “superbee”  limiter,  defined  as 

ip  2  —  max(0,  min(2i2, 1),  min(/2,  2)).  (4.46) 

This  limiter  introduces  a  large  artificial  compression  where  possible.  This  means  that 
discontinuities  in  the  solution  will  be  resolved  nicely,  however,  non-physical  sharp  layers 
can  appear. 

When  <P  =s  1,  rp+(R)  corresponds  to  another  Roe  limiter: 

rpi  =  max(0,min(R,l)).  (4.47) 

Note  that  0  <  ipi(R)  <  1.  This  means,  that  this  limiter  may  add  only  an  artificial 
viscosity  when  it  is  needed  to  damp  the  oscillations,  but  does  not  add  any  artificial 
compression.  The  solution  produced  by  S  scheme  using  this  limiter  will  contain  only 
layers  representing  the  physical  discontinuities.  However,  these  layers  will  not  be  as 
sharp  as  in  the  previous  case. 

Remark  4.6  Van  Albada  limiter  does  not  possess  the  property  (4.31),  therefore  if  used 
with  the  S  scheme,  a  non-monotone  solution  may  be  produced.  It’s  unique  property  is 

rl>VA  e  C°°.  (4.48) 

However,  it  is  not  needed  for  a  second  order  accuracy.  We  shall  also  show  that  this 
property  is  not  necessary  to  obtain  a  fast  convergence  of  the  multigrid  algorithm. 


4.2.5  Relaxation 


Denote 

F(u^)  =  h{fi+^  -  +  giJ+±  -  &j_$)  =  0.  (4.49) 

A  pointwise  relaxation  sweep  implies  updating  the  value  of  the  numerical  solution  utJ 
at  each  internal  grid  point.  This  can  be  done  performing  one  Newton  iteration  for  the 
nonlinear  equation  (4.49) 


(4.50) 
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4.3  2D  approach:  Inhomogeneous  equation 


To  implement  this  formula,  we  have  to  differentiate  with  respect  to  u,j  all  numerical 
fluxes  participating  in  the  balance  equation.  For  the  first  order  accurate  upstream  scheme 
these  derivatives  are  ^ 

(SUiJ+0mj 


=  a 
=  0 
=  i 
=  0. 


(4.51) 


For  the  first  order  N  schemer- 


t  (4.52) 

(Ai+*k 

(sVi)',. 

If  we  assume  differentiability  of  the  limiter  rf> ,  the  flux  derivatives  for  the  S  schemes  also 
appear  to  be  as  simple  as  for  (4.11)  (see  [20]).  They  are: 


=  ia 
=  0 

=  0. 


(/"**  A,  =  -  i(* 

=  (4.53) 

Newton’s  method  is  known  to  have  a  quadratic  convergence,  while  the  relaxation  process 
is  only  linearly  convergent,  even  locally.  Therefore,  relaxation  does  not  really  take  ad¬ 
vantage  of  the  fast  convergence  of  the  Newton  iterations.  This  is  the  reason  why  the  flux 
derivatives  can  be  calculated  approximately,  treating  the  limiter  function  as  a  quantity 
independent  of  u,  j.  Moreover,  the  flux  derivatives  of  the  S  scheme  can  be  substituted 
by  the  flux  derivatives  of  the  N  scheme.  The  numerical  experiments  confirm  that  this 
substitution  does  not  influence  the  performance  of  the  multigrid  solver. 

We  can  conclude  that  it  is  not  neccessary  for  limiters  to  be  C7  functions  for  the  fast 
multigrid  convergence  as  well  as  for  second  order  accuracy. 


4.3  2D  approach:  Inhomogeneous  equation 

The  S  scheme,  as  constructed  previously,  is  second  order  accurate  in  case  of  a  homoge¬ 
neous  equation.  We  «h»11  generalize  the  S  scheme  in  order  to  maintain  the  second  order 
accuracy  for  an  inhomogeneous  problem  (4.1). 


4.3.1  2D  scheme 


+*«-») 


Denote 


(4.54) 


Define 

~  ~  5(a(U*J  —  +  Kui-lj  ~  'J-l)  ~  K-Jj)  /.  __V 

97Dij-$  =  9Cij-$  ~  %W*ij  -  u.J-i)  +  a(uij_i  -  Ui.u^x)  -  hs{J_i) 

Lemma  4.4  Scheme  (4.55)  is  second  order  accurate. 

Proof:  This  lemma  can  be  proved  in  a  similar  way  to  the  corresponding  homogeneous 
case.  We  want  to  show  that 

r°»hJ-rD<-u  -««.)+. o(fc»)  ■ 

s’Vi  -  »JVi  =  MK)  +  0(*3)-  ( 5  ) 

Note,  that 

f'D^J-f2Di-U  =  fC^~fCi-U- 

~  «.j)  +  6(uij  -  Uij.O  -  K+*j)  + 

|(«(u»j  -  tti-u)  +  6(tn_ij  -  Ui-w-i)  -  hs^j)  (4.57) 

Since  the  central  scheme  (4.3)  is  second  order  accurate,  we  can  conclude  that 

f2D_  .  _  f2D  _ 

A(a«s)  —  iA2(eti,  +  An*  —  s),  +  0(A3),  (4.58) 

or  taking  the  entire  equation  into  account 

-  /“.-*„•= 

=  *(«“.)  +  0(fc3).  (4.59) 

Similarly 

92Dij+l  -  =  *(K)  “  5**(a«»  +  K  -  ■»)*  +  ^(h3)  = 

A(K)  +  0(*8).  (4-60) 


(4.57) 


(4.58) 


(4.59) 


(4.60) 


4.3.2  N  scheme 

In  case  one  is  interested  to  obtain  first  order  accuracy  the  previously  defined  N  scheme 
can  be  used.  However,  the  purpose  of  it  here  is  to  serve  as  an  intermediate  step  towards 
the  construction  of  the  second  order  accurate  S  scheme.  Therefore,  we  shall  modify  it. 
For  our  representative  case  a  <  b,  putting 

0  =  min(l,^)  =  ^  <  1,  (4.61) 


the  N  scheme  can  be  defined  by: 

/V jj  =  ~  “  “*-»•>)  +  PW*- U  -  ««- W-i)  - 

9*ij-}  =  sfij-j  -  iW&j  ~  «u-i)  +  -  «i-w-0  -  Km)- 


(4.62) 
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4.3.3  S  scheme 


The  S  scheme  is  defined  by 
f*i-y  =  /V y  - 

*Vj  =«Vi-  -  -  •  • 

-  —  .  .  •  i“ 

whore 

i?.  ' 

*“*J  b{ui-U  ~  ~  hsi.^j 

Lemma  4.5  If  V>  —  V»(i2)  €  C*  and 

*(1)  «  1, 

then  the  S  scheme  is  second  order  accurate. 

Proof:  In  order  to  prove  this  lemma  it  is  sufficient  to  show  that 


(4.63) 
.  (4.64) 


(4.65) 


=  +  (4.66) 

Observe  that 

/??*„•  -  f3,-y  =  5(1  -  WRi-iM1  ~  MKv-u  ~  “i-w-0  -  H-*.)  (4-67) 

Taking  into  account: 


1  -  *(*-*  j)  =  -  1)  +  0(h3)  (4.68) 

Then 

/2P  rS  . 

-^k(i)(«i-4j  -  i)(i 

Therefore 

/»*  -  /VjJ  =  Cj  J  -  /“(-J.,  +  0(h\  (4.69) 

which  together  with  the  fact  that  2D  scheme  is  second  order  accurate  proves  (4.66). 


-  + o(h3) = 

-  u— u-i) + -  “i-i j-t) + + °(*3)- 


Remark  4.7  In  order  to  obtain  a  second  order  accurate  solution  to  the  inhomogeneous 
problem  by  the  FMG  algorithm  it  is  neccessary  to  use  this  scheme  on  the  currently 
finest  grid  only.  The  scheme  constructed  for  homogeneous  problems  can  be  employed 
on  coarse  grids.  The  coarse  grid  correction  obtained  this  way  will  be  only  first  order 
approximation  to  the  needed  correction  on  the  finest  grid.  However,  it  is  satisfactory, 
because  the  needed  correction  is  only  0(h7)  large. 


Chapter  5 


Discretization:  General  case 


There  are  several  slightly  different  ways  to  extend  the  genuinely  2D  approach  for  the 
case  of  a  nonlinear  conservation  law.  We  shall  present  here  one  of  the  simplest.  We  shall 
first  extend  the  construction  of  the  difference  scheme  and  prove  its  monotonicity  for  a 
homogeneous  case.  Then  we  shall  treat  an  inhomogeneous  case,  demonstrating,  that  a 
second  order  accuracy  can  be  obtained  by  this  approach.  A  second  order  accuracy  of  the 
schemes  constructed  for  a  homogeneous  case  will  follows  directly  from  the  more  general 
result. 


5.1  Homogeneors  equation 

Consider  a  2D  homogeneous  conservation  law 


(/(«)).  +  (sM),  =  «■ 

(5.1) 

Denote 

•4  <5 

II  II 

"3'^S' 

(5.2) 

Then  (5.1) 

can  be  written  in  quasilinear  form 

a(u)us  -f  6(u)uy  =  0. 

(5.3) 

5.1.1  Central  and  upstream  schemes 

A  central  unstable  second  order  accurate  difference  approximation  for  (5.1)  is  determined 
by 


=  j(/iJ  +  /*-  w) 

9°  =  jG/ij  +  9iJ- 1)’ 

In  order  to  stabilize  (5.4),  we  add  artificial  viscosity  term 


S f  i— ala*— f  ***'"W) 


(5.4) 


(5.5) 
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where 


and 


fD<-y 


=  fci-U  -  -  u«-u)  +  4>t\j  +  Ky) 

=  9‘iJ-h  ~  +  7+-_*  +  \M), 


(5.14) 


7^-1  =  »gn(JBi+iJ_i)A-t+,  (5xu)i+^_^. 


The  second  order  accuracy  of  this  scheme  will  be  a  direct  corollary  from  that  of 
the  more  general  scheme,  approximating  an  inhomogeneous  equation,  which  will  be  con¬ 
structed  in  Sec.5.2.1.  However,  our  main  concern  now  is  monotonidty.  The  scheme  given 
by(5.14)  is  not  monotonic. 


5.1.3  First  order  N  scheme 

Denote 


Modify  the  previoiisly  defined  upstream  scheme 

f"i-y  =  fei-y  -  -  «.-w) 

-  Uij-i)' 


(5.17) 


(5.18) 


Note  that  the  downstream  grid  points  may  also  participate  with  positive  coefficients 
in  the  difference  equation  defined  by  (5.18).  However,  they  will  be  only  0(h)  large 
comparing  with  coefficients  in  the  upstream  grid  points.  Therefore,  relaxation  sweeps  in 
the  downstream  direction  will  still  be  very  efficient  way  to  solve  the  difference  equations. 
Denote 

a+i-y-i  =  min(l, 

Q~i-y+k  =  min(l, 


ifbiizii) 

A+H  j-i 

3-i.jil1 

A* .  i  . .  i 


) 


(5.19) 


and 


P  i+ii-k 


—y-i 


—  min(l,  ■ 


). 


We  shall  call  -“N  scheme”  the  first  order  scheme  defined  by 


(5.20) 


5.1  Homogeneous  equation 


=  9U  ~  t  7 

(5.21) 

where 

. 

- 

*£i*  =  j  -  . 

CjJ  -  *  ;  - 

■  (5.22) 

and 

Vi  = 

(5.23) 

It  is  easy  to  see  that  in  case  of  a  linear  constant  coefficient  equation  this  scheme 
coincides  with  the  N  scheme  defined  in  the  previous  chapter. 

Theorem.  5.1  The  N  scheme  defined  by  (5.21)  is  of  positive  type. 


Proof:  When  (5.21)  is  substituted  into  the  balance  equation,  the  resulting  relation  can 
be  written  in  the  following  form: 

X)  CkAUi+k’i+l  ~  u*j)  “  (5.24J 

W+M/o 

where  Jfc,  /  =  —1, 0, 1.  We  want  to  show  that  all  the  coefficients  C*  j  >  0  for  |£|  +  |/|  /  0. 
In  order  to  show  this,  we  have  to  consider  all  the  possible  situations  for  each  gridpoint 
participating  in  the  equation  and  to  verify  the  non-negativity  of  the  corresponding  coef¬ 
ficient. 

Consider  the  diagonal  point  i  —  l,j  —  1 .  The  corresponding  coefficient  will  be  deter¬ 
mined  by  the  quantities  calculated  for  the  grid  square  i  —  5,  j  — 


CV.  = 


(5.25) 


or 


C:i,_1=min(A^..x,B+,j..p>0 


(5.26) 


Non-negativity  of  other  coefficients  corresponding  to  the  diagonal  points  can  be  shown 
in  a  similar  way. 


Consider  the  gridpoint  i,j  —  1. 


<*-. = £(Vi + Vi  -  «?-ij-i  - 

(5.27) 

where 

.  i 

(5.28) 
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Recalling 

it  is  easy  to  see  that 

C0°.,i  >  0. 

Non-negativity  of  other  coefficients  can  be  verified  in  a  similar  way. 


(5.29) 


(5.30) 


5.1.4  SI  and  S2  schemes 


We  want  to  add  a  second  order  correction  to  the  N  scheme,  maintaining  its  monotonidty. 
We  shall  do  it  in  two  different  ways  (Si  and  S2).  Two  difference  schemes  will  be  con¬ 
structed.  One  of  them  (the  Si  scheme)  is  slightly  simpler  and  has  the  smaller  truncation 
error,  however  it  may  admit  non-physical  discontinuities  and  is  proven  to  be  monotonic 
only  when  used  with  non-compressive  limiters  (like  Roe’s  ij> i  limiter).  Another  scheme 
(S2)  introduces  larger  numerical  viscosity,  but  rejects  non-physical  discontinuities  and  is 
proven  to  be  monotonic  when  used  with  compressive  limiters  as  well.  Define 


(5.31) 


=  (1  -  °”  »+$  )s*8uC®»'+$.i- $  WiQi+l  j-J  £  (£*“)„•+$ 


where 


(5.32) 


(5.33) 


Then  the  Si  scheme  can  be  defined  by 

f*t-U  =  /V  H  -  +  ♦“'H) 

»“u-i  = 

We  can  formulate  the  following  monotonidty  result 


(5.34) 


(5.35) 


Theorem  5.2  The  SI  scheme  defined  by  (5.35)  is  monotonic  if 


o  S  <  1,  V’(iJ)  <  1. 


(5.36) 


5.1  Homogeneous  equation 
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Proof:  Using  the  following  identities 


«2zH-k] 

Q-*-M 


(5.37) 


it  can  be  verified  that  it  is  always  possible  to  rewrite  the.Sl  scheme  in  the  form 


2  &kj(ui+kj+i  ~  Uij)  —  0,  (5.38) 

W+HI#0 


where  fc,/  =  —1,0, 1  and,  provided  (5.36)  holds  Cj$  >  0  for  |fc|  + |/|  ^  0. 

The  complete  proof  is  straightforward.  We  shall  only  illustrate  the  importance  of  the 
restriction  (5.36)  on  two  typical  situations,  considering  grid  point  (:,  j  —  1). 


1.  Suppose 


Ai-iJ-i  =  >  0 

(5.39) 

This  means  that 

(*v =  «i-xj  ~ 

(5.40) 

Suppose  also 

B{+jJ-}  = 

(5.41) 

and  that  A^j 

is  a  small  positive  number.  Then 

=  Ui+ij-i  ~  «y-i 
(^VU)i+$j-$  “  ViJ  ~  ViJ- X* 

(5.42) 

This  yields 

ctx  =  CS.-X  -  (1  - 

(5.43) 

or 

^0,-1  —  **■  ~  ~  A«+i<i-§ 

-(4+iJ-i  -  A+iJ-Om+U-i))  (5.44) 


The  last  expression  will  be  non-negative  if  (5.36)  holds,  but  may  attain  negative 
values  if  tff(R)  >  1. 


2. 


Suppose  that  >  0  and  that  <  0  has  a  small  absolute  value.  This 

means  that 


=  wy-i  ”  «i-W- 1 
=  Wj  ~  «W- 1- 


(5.45) 
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(5.47) 


Suppose  also 

(5‘46) 

and  is  a  small  positive  number.  Then 

■  .  ...  (5.47) 

Note  that  this  case  corresponds  to  a  rarefaction  wave.  Here  we  obtain 

-(1  -  (5.48) 

Since  and  are  small,  the  last  expression  may  attain  negative  values 

if  >  1.  However,  it  will  be  non-negative,  if  (5.36)  holds.  If  we  rewrite  the  Si 
scheme  corrections  using  the  identities  (5.37),  it  will  create  negative  coefficients  at 
diagonal  grid  points  (*  —  1  ,j  —  1)  and  ( i  +  1,  j  —  1). 

Remark  5.1  The  inequality  (5.36)  means  that  no  artificial  compression  can  be  added 
to  the  schemes. 

Remark  5.2  The  SI  scheme  (as  well  as  the  N  scheme)  will  perfectly  resolve  contact 
discontinuities  and  shock  waves  which  align  either  with  grid  lines  or  with  grid  diagonals. 
However,  non-physical  discontinuities  (corresponding  to  rarefaction  waves)  can  also  be 
admitted  in  these  cases.  This  is  because  of  the  vanishing  cross-stream  truncation  error. 

We  shall  now  modify  the  Si  scheme  in  order  to  enable  it  to  reject  non-physical 
discontinuities  and  to  allow  the  addition  of  artificial  compression  when  needed.  Denote 


(«,*), •  (5.49) 

where  n  and  v  are  the  same  as  in  definitions  of  and  (5.9).  Denote  also 

1  =  ««  -  «i-W  ~  «.J-i  +  «i-W- 1  (5-5°) 


=  2 -  Ui-u)  -  (e*A)i_iii_i(^u)i_iJ_| 

ij  =  2(^B),_iJ_i(uij  -  Uf-u)  +  (foA)i_iti+i(^u)i.iJ+i 


T+U-i  = 


2(feA)f_i>J_$(ti,v  -  Uij-0  -  (eM-ij-tiQuh-jj-t 
2(PrA){+jj_x(u,j  -  ufJ_i)  +  (evB){+ij_i(Slyu)i+ij_i 


(5.51) 


(5.52) 


Define 


(5.53) 


Theorem  5.3  The  S2  scheme  is  monotonic  if 

0  <  <  2,  1>(R)  <  2.  (5.54) 

Jti 

Proof:  Once  again  we  claim that  it  is  possible  to  rewrite  the  S  scheme  using  the  identities 
(5.37),  wheremeccessary,  in  the  following  form  ~ 

-  -  ^2  -:  Ckj(ui+kj+i  —  uijh=  0,  ^  .  (5.55) 

-  l*l+M*o 

where  k,l  =  —1,0, 1  and  CfJ  >  0  for  \k\  +  |/|  ^  0. 

We  shall  decsribe  in  details  only  the  same  two  possible  situations,  as  in  the  previous 
theorem,  considering  a  grid  point  (i,j  —  1): 


1.  Suppose 

>  ®.  (5  56) 

J  =  Bi-iJ-i  (557) 

and  a  small  positive  number. 

Then 

^,  =  C|Ll  +  ((fcA).-.lj-J+(foA),+jJ..})  (5.58) 

Recalling  the  expression  for  Co,-i  and  taking  into  account  the  following 

((exA)i_^_i  +  (e,i4)i+jj_j)  >  -  A+ij-ili  (5-59) 

we  can  conclude  that  C$t i1  is  nonnegative  if  (5.54)  holds. 

2.  Suppose 

5<+id-i  =  (5-6°) 

and  is  a  small  positive  number. 

Suppose  also  thatl  >  0  and  <  0  has  a  small  absolute  valuer. 

Note  that,  since  and  A+Jj-J  have  opposite  signs,  then  either 

(f,A)HiM  >  (5.61) 


or 

(foA)i+jj_j  >  (5-62) 

Suppose  for  simplicity,  that  the  first  is  correct.  Then  the  SI  scheme  correction 
corresponding  to  »  •  »  and  B{  i  •  i  can  be  substituted  by  using  (5.37)  and  we 
shall  obtain 


Cff..,  =  (ft*),-*.*  + 


(5.63) 
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+  hj-i 

'  :  '  +(1  . 

It  is  easy  to  see  that  both  C£f  and  C^lx  will  be  nonnegative  if  (5.54)  holds. 

Remark  5.3  The  S2  scheme  which  was  created  from  the  Si  scheme  by  adding  the  e- 
correction  allows  us  to  use  limiters  which  may  introduce  artificial  compression. 

Remark  5.4  Another  important  property  of  the  e-correction  is,  that  it  will  introduce 
additional  cross-stream  viscosity  in  case  of  a  rarefaction  wave  (diverging  diaracteristic 
field).  Therefore,  no  non-physical  discontinuities  will  be  admitted. 

5.2  Inhomogeneous  equation 

Consider  an  inhomogeneous  conservation  law 


(/(«)),  +  (s(u))v  =  a, 


(5.64) 


where  s  =:  s(x,  y). 


5.2.1  2D  scheme 

The  same  central  and  upstream  approximations  can  be  used  in  this  case  and  they  will 
be  second  and  first  order  accurate  respectively.  However,  the  genuinely  2D  scheme  has 
to  be  modified,  in  order  to  maintain  the  second  order  accuracy. 


Denote 


si-$j  =  +  *•- la) 

3ij-\  ~  5*Vi-0 


Define 


.+  .  - 


(5.65) 


(5.66) 


(5.67) 


S.2  Tn homo geneoua  equation. 


*-*  AU,-i  -  Ar+i,-} 


at.  i  = 


(5.68) 


Again 


-  -  • 

rVi,  =  /*Ma  -  iW-j, + Cj.) 

S^J-i  =  +  !«-))■ 

^  =  riga(/t,.i^i)(B+.-^-)(«»“).-)J-i  +  h3UJ 

$-_y  =  sign(Xi.jiJ+j)(B"i.jJ+j(J,“)j-)j+)  +  ^i-}  j) 

7+  i  =  +  K.-P 

7"._,  ss  sigu(Bt>^_j)(A”i+jj.„j(6«u)i+^_j  + 


(5.69) 


(5.70) 


(5.71) 


(5.72) 


Theorem 

accurate. 


5.4  He  2D  scheme  defined  by  (5.70)  wiiA  (5.71)  end  (5.72)  is  aaonsi  »<*r 


Proof:  Assume  that 


3ign(A_ij-i)  =  sWA-iJ+i)  =  sign(At+iii_p  -  sign(A,-+*j+*) 

sign(B,_^_p  =  sign  (Bi-jj+j)  =  sign(Bi+^_p  =  sign(Pi+|j+|). 


Then 


as  6aaign(a(u))(6(u)u,  —  a)*  +  0(5*). 


Assuming  that 


we  also  obtain 


•W4-W+J  -  A-lj+j)  =  ®P>(^Ja+J  ~ 

1 

-  i(*Jngn(«(ii))(a(ti)u.),  +  **•« o(«, )(«,«.).)  + 0(5S) 

it 


(5.73) 


(5.74) 


(5.75) 


(5.76) 


Since  the  central  scheme  is  second  order  accurate 


=  (/(«)).  +  0(k°) 

+|(jai+^i(u»+w  -  -  (tf+ij + tr+iJ) 

Taking  into  account  the  previous  three  relations,  we  conclude 


(5.77) 


tlD  r2D 

■'*-  2  J 

=  (/(«))*  +  h2sign(<x(u))(a(u)ur  +  b(u)uy  -  s)x  +  0(h 3)  =  0(h3)  (5.78) 

Similarly 

2D  2D 

=  (y(u))x  +  h2sign(6(u))(a(u)t«x  4-  &(u)u„  -  s)„  +  0{hz)  =  0(h3)  (5.79) 

If  one  of  the  assumptions  (5.73)  or  (5.75)  does  not  hold,  the  difference  equation  in  this 
computational  cell  will  approximate  the  differential  one  only  up  to  0(h).  However,  this 
can  happen  in  computational  cells  which  cover  only  0(h)  part  of  the  domain.  Therefore, 
the  scheme  will  be  still  second  order  accurate. 


5.2.2  N  scheme 


If  we  are  interested  in  obtaining  first  order  accuracy  only,  we  can  use  the  same  N  scheme, 
as  for  the  homogeneous  case.  However,  if  the  N  scheme  is  an  intermediate  step  towards 
constructing  higher  order  schemes,  it  also  has  to  be  modified. 

Once  again 


where 


ft-iJ 

9° 

=  /Vi,  - 

=  ~  a(<7°+»,7-f  7°  ij-j) 

-  (5.80) 

(5.81) 

(5.82) 

5.2.3  S  scheme 


Define 

At  =- 

B.*-iJ+i(^4|)*-^+l +  A3‘*-iu 

Ot. 

ij_*  +  **y-J 

Q_  _ 

Then,  the  Si  scheme  can  be  given  by 


(5.83) 

(5.84) 

(5.85) 

(5.86) 


=  fi-ij  -  K^Vjj  + 

«*!«-»  +  «-}) 

with 

^+._^  = 

(1-  P~  +  4<-ij) 


(5.87) 


(5.88) 


751+«a-i  = 

(X-  jJsign^.^XA.j  +  ««-*)  (5  g9) 

"751  fj-J  = 

(1-  a-i+ij_i)^(Q-^i)sign(Bi+^_i)(A+4j_4(^«),+4j-4  +««-§)• 

Theorem  5.5  If  r}>  =  ^(i2)  €  C3  and 

*(1)  =  1,  (5.90) 

then  the  SI  scheme  is  second  order  accurate. 

Proof:  In  addition  to  (5.75)  and  (5.73),  assume  also  that 

>  0.  (5.91) 
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This  means  that 


Btu-i  =  >  0 

=0’ 

Cur?  . 

'  € -u-»  . 


‘7-y =  0 

<U  =  sH-r 


(5.92) 


(5.93) 


(5.94) 


Then 


(2D  fs\  _ 

A-iJ  “  A-iJ  = 


-50  - 

Taking  into  account 

1  -  «A)  =  tf,(l)(A  - 1)  +  0(fc3), 


(5.95) 


(5.96) 


we  have 

(2D  fs\  _ 

A-iJ  A-Jd 

+  0(h9)  (5.97) 

Recalling  the  definition  of  A+  we  can  rewrite  (5.97)  as 

fXD  /SI  _ 

A-iJ  “  A-iJ  ~ 

•W^-U-jX1  -  f>tu)  +  0(*3). 

Then 

/"u  ~  J?u  -  "  ^?U + (5-98> 

If  one  of  the  assumptions  (5.73),  (5.75)  or  (5.91)  does  not  hold,  the  difference  equation 
at  that  computational  cell  will  approximate  the  differential  one  only  upto  0(h).  How¬ 
ever,  this  can  happen  in  computational  cells  which  cover  only  0(h)  part  of  the  domain. 
Therefore,  the  scheme  will  be  still  second  order  accurate. 


Remark  5.5  The  S2  scheme  created  by  adding  the  e-correction  to  the  SI  scheme  will 
also  be  second  order  accurate.  This  is  because  the  e-correction  is  0(fc3). 


Chapter  6 

Numerical  experiments 


All  the  numerical  experiments  reported  here  deal  with  numerical  solutions  of  the  differ¬ 
ential  equation 

-  (0+)Au  +  (/(u))x  +  (g{u))y  =  3,  (6.1) 

where  0+  means  an  infinitesimally  small  positive  number. 

We  choose  for  our  domain  the  rectangle  : 

a  =  {(x,y):0<x<3,0<y<2}  (6.2) 

and  we  set  Dirichlet  boundary  conditions  around  its  boundary. 

We  used  five  qrids  (levels).  The  meshsize  of  grid  k  is  h*  =  21-*,  hence  it  has  (3  x 
2*_l  —  1)  x  (2*  —  1)  interior  grid  points  and  5x2*  boundary  points  (1  <  k  <  5). 

Note  that  in  all  the  cases  reported  below  limit  solutions  can  be  obtained  just  by 
few  downstream  relaxation  sweeps  on  the  finest  grid.  This  is  because  all  the  difference 
schemes  used  are  upstream  (or  “almost”  upstream).  However,  Eq.(6.1)  is  just  a  model 
problem  for  more  complicated  systems  which  cannot  be  solved  this  way.  In  case  of  sub¬ 
sonic  flow  for  instance  the  equations  contain  an  elliptic  component  as  well  as  hyperbolic. 
In  case  of  supersonic  flow  there  exist  several  families  of  characteristics.  Therefore,  the 
main  purpose  of  the  experiments  reported  here  is  to  demonstrate  that  the  developed 
discretization  schemes  provide  the  possibility  to  obtain  second  order  accurate  solution 
(in  smooth  regions  and  also  in  terms  of  discontinuity  location)  by  means  of  a  certain 
multigrid  algoritm,  employing  a  direction-free  relaxation. 

The  algorithm  used  is  of  type  FMG(N ,  JVm,  (7,  Af)  (see  Chap.3.),  where  M  =  5,N  = 
1,2  ,Nm  as  2,6,  and  C  =  W( 2,1).  The  Full- Weighted  residual  transfer  is  used.  In 
case  the  N,  S,  SI  or  S2  schemes  are  used,  the  Red-Black  relaxation  without  storing  of 
intermediate  values  is  not  direction-free  anymore.  Therefore,  we  use  “4-colour”  ordering. 
The  usual  bilinear  correction  interpolation  and  bicubic  FMG  interpolation  is  employed. 

The  precise  formulas  for  numerical  flux  derivatives  were  used  in  all  the  experiments 
reported  here.  However,  there  will  not  be  any  significant  difference  in  the  performance 
of  the  algorithms  if  the  N  scheme  numerical  flux  derivative  formulas  will  be  used  for  the 
SySl  and  S2  schemes.  This  has  been  observed  in  a  comparison  of  the  solution  error  and 
residual  behaviour  in  both  cases. 

We  shall  compare  the  solutions  obtained  by  different  discretizations  and  discuss  the 
choice  of  the  difference  scheme  for  a  particular  problem. 
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6.1  Linear  equation 

the  quality  of  obtained 
-  (6.3) 

6.1.1  Smooth  solution'  - 

First  we  study  how  effective  different  algorithms  are  in  case  of  smooth  solutions,  avoiding 
any  influence  of  discontinuities.  For  this  purpose  we  provide  equation  (6.3)  with  such 
boundary  conditions  that  there  will  be  no  boundary  layers.  Since  all  the  difference 
schemes  we  experiment  with  here  are  upstream,  we  do  not  expect  the  appearance  of 
even  a  numerical  boundary  layer. 


We  shall  first  examine  the  performance  of  the  algorithms  and 
numerical  solutions  in  case  of  the  linear  equation 

—  (D+)Au  +  an*  +  bu^  as  s. 


Homogeneous  case 

Consider  the  following  version  of  equation  (6.3) 

—  (0+) Au  +  .5u*  4-  «y  »  0, 

(6.4) 

with  boundary  conditions  given  by  ;  *-■  Y 

«  =  sin(y  —  2x). 

(6.6) 

It  is  easy  to  see  that  (6.5)  is  also  the  exact  solution  of  (6.4).  This  solution  is  constant 
along  the  characteristics  and  varies  in  other  directions.  We  want  first  to  compare  per¬ 
formance  of  the  upstream,  N,  2D  and  S  schemes  on  this  model  problem.  The  S  scheme 
is  employed  here  in  three  versions:  with  Van  Leer’s  limiter  and  two  Roe’s  limiters.  The 
experiments  with  this  model  problem  are  presented  in  Table  6.1. 

Each  column  of  Table  6.1  starting  from  the  second  presents  the  history  of  the  Li 
norm  of  the  solution  error  (the  difference  between  the  numerical  and  the  differential 
solution).  2 FMG  algorithm  is  employed  (i.e.,  N  =  2).  The  first  column  indicates  the 
stage  of  the  multigrid  algorithm  the  error  is  displayed  at:  the  first  number  stands  for 
the  currently  finest  level  and  the  number  in  parentheses  says  how  many  multigrid  cycles 
have  already  been  performed  on  this  level  (zero  means  that  the  error  is  displayed  just 
after  the  bicubic  interpolation  to  this  level).  Columns  1,2  correspond  to  the  cases  where 
the  upstream  and  N  schemes  are  used  respectively.  Both  these  schemes  demonstrate  first 
order  convergence.  However,  the  solution  error  in  case  of  the  N  scheme  is  almost  three 
times  smaller.  The  IF  MG  algorithm  will  produce  results  of  the  same  quality  in  both 
these  cases  (This  can  be  seen  ,  e.g.,  by  comparing  the  results  of  row  5(1)  with  those  of 
5(2)  and  5(6);  the  latter  practically  show  the  discretization  error).  Column  3  presents  the 
experiments  with  the  2D  scheme  and  Columns  4,5  and  6  correspond  to  the  experiments 
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Difference 

scheme 

Upstr 

N 

2D 

S 

S 

S 

Limiter 

- 

- 

- 

1>vl 

fa 

fa 

2(5) 

.241 

.119 

.0429 

.0497 

.0588 

.0471 

-  3(2)  ' 

:  .154 

.0632 

.00971 

.0145 

6232 

.-6101 

4(2) 

.0870 

.  .0326 

.00251 

.00398 

60688 

.00638 

.  5(0) 

.0846 

.0317 

.00243 

.00376 

.00656 

.00597 

5(1) 

.0475 

-.0161 

.00054 

.60150- 

60266 

60275 

5(2) 

.0467 

.0165 

.00065 

60113 

.00213 

.00235 

5(6) 

.0468 

.0166 

.00060 

.00097 

.00194 

.00219 

1 

r  2 

3 

4 

5 

6 

Table  6.1:  Linear  homogeneous  problem;  slowly  varying  solution. 

with  Van  Leer’s  limiter  and  Roe’s  fa  and  fa  limiters  respectively.  The  second  order 
convergence  can  be  observed  in  all  these  cases. 

Table  6.2  presents  the  experiments  with  the  same  equation  (6.4)  but  with  different, 
more  oscillatory  boundary  conditions: 

t*  =  sin(5(y  —  2x))  (6.6) 

and  has  the  same  structure  as  Table  6.1. 


Difference 

scheme 

Upstr 

N 

2D 

S 

S 

S 

Limiter 

- 

- 

- 

fan. 

fa 

fa 

2(5) 

.721 

.759 

.972 

.779 

.771 

.78 8 

3(2) 

.644 

.577 

.860 

.557 

.561 

.553 

4(2) 

.584 

.434 

.382 

.267 

.317 

.202 

5(0) 

.569 

.429 

.358 

.274 

.322 

.208 

5(1) 

.515 

.327 

.131 

.0994 

.155 

.0806 

5(2) 

.496 

.298 

.123 

.0728 

.121 

.0705 

5(6) 

.490 

.294 

.0835 

.0675 

.113 

.0542 

1 

1 

2 

3 

4 

5 

6  1 

Table  6.2:  Linear  homogeneous  problem;  oscillatory  solution. 


The  upstream  and  N  schemes  seem  to  give  a  very  poor  approximation  to  such  a 
solution,  however  the  2D  and  S  schemes  with  all  the  limiters  start  to  demonstrate  a 
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second  order  convergence  on  the  finer  levels.  However,  about  3  cycles  are  needed  on  the 
finest  level  to  achieve  a  second  order  accurate  solution.  Performing  more  than  2  cycles  on 
each  coarser  level  does  net  change  the  situation,  because  the  coarse  grid  solution  provides 
a  poor  approximation  for  the  finer  grid  one  (for  detailed  analysis  of  this  situation  see 
Sec.2.3  in  [4]). 


Inhomogeneous  ease 

Consider  the  following  inhomogeneous  equation 

—  (0+)  Au  +  .5U*  +  tty  =  -s, 

with  boundary  conditions 


Choose  the  right-hand-side  to  be 


u  =  sin(x  +  y). 


s  =  1.5cos(x  +  y). 


Then  (6.8)  will  be  the  exact  solution  of  (6.7). 

Table  6.3  has  also  the  same  structure  as  Table  6.1. 


Difference  II  N  N 

scheme  I  Upstr  horn,  inhom.  2D 


Limiter 


s 

S 

s 

fa 

fa 

.165 

.0844 


.107 

.0277 


.0414  .0717  .00957  .00765  .00596  .00562 


.0347 

.00393 

2 

3 

.00136  .00231 


Table  6.3:  Linear  inhomogeneous  problem. 

Once  again  the  algorithm  employing  the  upstream  scheme  demonstrates  a  first  order 
convergence,  as  resulting  from  Column  1.  Columns  2  and  3  present  the  experiments 
with  two  versions  of  the  N  scheme  -  without  weighting  of  the  right-hand- side  (we  shall 
call  it  "homogeneous”)  and  inhomogeneous  (which  was  constructed  as  an  intermediate 
stage  towards  the  second  order  accurate  S  scheme,  see  Sec.4.3).  The  homogeneous  N 
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scheme  demonstrates  first  order  convergence  and  its  solution  error  is  even  larger  than 
that  of  the  upstream  scheme.  However,  the  inhomogeneous  N  scheme  leads  to  the  much 
smaller  error  and  even  seem  to  demonstrate  the  higher  order  convergence  on  coarse  grids. 
This  is  because  the  streamwise  component  is  second  order  small  in  this  case  (due  to  the 
weighting  of  the  right-hand-side),  but  it  can  dominate  the  cross-stream  first  order  error 
on  some  coarse  grids.  The  algorithms  based  upon  the  2D  and  S  schemes  (modified  for  the 
inhomogeneous  case)  using  anyone  of'the  three  limiters  produce  second  order- accurate 
solutions.  In  the  cases  of  the  2D  and  S  schemes  with  Van  Leer’s  and  Roe’s  ^  limiters 
this  is  already  clearly  achieved  by  IF  MG  algorithm.  When  Roe’s  ^  limiter  is  employed 
the  convergence  may  be  just  a  little  slower. 

6.1.2  Resolution  of  contact  discontinuities 

We  shall  examine  the  performances  of  these  algorithms  in  case  of  contact  discontinuity. 
Consider  Eq.(6.4)  with  boundary  conditions  given  by 

u  =  H( 0.5y  -  x  +  1),  (6.10) 

where  H  is  the  Heaviside  function:  H(x)  =  0  for  x  <  0,  H{x)  =  1  for  x  >  0. 

It  is  easy  to  see  that  (6.10)  is  also  the  exact  solution  of  (6.4)  under  these  boundary 
conditions,  and  it  contains  a  jump  discontinuity  along  the  line  y  =  2x  —  2. 

Figures  6.1  (a— /)  present  the  numerical  solutions  to  this  problem  obtained  by  2FMG 
algorithms  employing  different  discretization  schemes.  The  results  of  the  IF  MG  algo¬ 
rithm  can  hardly  be  distinguished  from  those  of  2 FMG  algorithm  (except  the  case  of  the 
S  scheme  with  Roe’s  fa  limiter),  therefore  we  omit  them.  Figures  6.2(o  —  /)  correspond 
to  the  same  numerical  experiments,  but  present  a  plot  of  the  solution  along  the  gridline 
y  =  1.75  for  1  <  x  <  3.  The  line  with  “cross”  points'  represents  solutions  obtained  on 
level  4,  and  the  line  with  “diamond”  points  represents  the  level  5  solutions  -  the  same 
as  in  Figure  6.1. 

Figures  6.1a,  6.2 a  and  6.16,  6.26  correspond  to  the  upstream  and  N  schemes  respec¬ 
tively.  The  contact  discontinuity  is  resolved  better  in  case  of  the  N  scheme.  However, 
this  result  is  still  unsatisfactory  because  the  width  of  the  transition  layer  decreases  only 
by  factor  y/2  when  the  grid  becomes  twice  finer,  i.e.,  the  number  of  gridpointa  in  the 
transition  layer  increases  by  roughly  a  factor  of  V2,  as  can  be  observed  in  Figure  6.26. 
Figures  6.1c  and  6.2c  correspond  to  the  2D  scheme.  The  spurious  oscillations  can  be 
observed  in  the  solution. 

Figures  6.1(d  —  /)  and  6.2(d  —  /)  correspond  to  the  S  scheme  using  Van  Leer’s 
and  Roe’s  ^  and  tfa  limiters  respectively.  The  transition  layer  is  slightly  wider  in 
case  of  Roe’s  than  Van  Leer’s  limiter.  The  discontinuity  profile  produced  by  the 
S  scheme  with  Roe’s  limiter  did  not  converge  enough.  This  is  because  the  large 
amount  of  artificial  compression  is  introduced  by  such  a  limiter  in  the  neighborhood  of  the 
discontinuity,  which  badly  affects  the  smoothing  properties  of  the  scheme.  The  solution 
can  be  improved  by  local  relaxation  sweeps  in  the  neighborhood  of  the  discontinuity. 


Figure  6.1:  Contact  discontinuity;  (a)  -  Upstream  scheme,  (b)  -  N  scheme. 
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Figure  0.2;  Contact  discontinuity;  gridline  y  =  1.25  for  1  <  x  <  3,  x  -  level  4,  o  - 
level  5:  (a)  -  Upstream  scheme,  (b)  *  N  scheme,  (c)  -  2D  scheme,  (d)  -  S  scheme  with 
Van  Leer’s  limiter. 


Figure  6.2:  continued;  S  scheme  with:  (e)  -  Hoe’s  rfo.  limiter,  (f)  -  Roe’s  tf> j  limiter,  (g) 
-  Roe’s  2  limiter,  limit  solution. 
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Figure  6.3:  Contact  discontinuity;  (a)  -  Upstream  scheme,  (b)  -  N  scheme;  S  scheme 
with:  (d)  -  Van  Leer’s  limiter,  (e)  -  Roe’s  fa  limiter,  (g)  -  Roe’s  fa  limiter,  the  limit 
solution. 

Figures  6.1$r  and  6.2 g  present  the  numerical  solution  to  the  same  problem  obtained 
by  few  relaxation  sweeps  in  downstream  direction  on  the  finest  grid  using  the  S  scheme 
with  Roe’s  fa  limiter.  The  transition  layer  created  is  very  sharp.  It  is  possible  to  obtain 
the  same  result  by  one  downstream  relaxation  sweep,  just  performing  several  Newton 
iterations  at  each  grid  point.  Note  that  downstream  relaxation  sweeps  can  be  performed 
locally  in  the  neighborhood  of  the  discontinuity  together  with  direction  free  relaxation 
all  over  the  domain. 

Figure  6.3  presents  together  the  same  (level  5)  numerical  solutions  which  appear  on 
Figures  6.2a,  6,  d,  e,  g. 

The  sharpest  discontinuity  profile  here  corresponds  to  the  limit  solution  of  the  S 
scheme  with  Roe’s  fa  limiter.  The  S  scheme  with  Van  Leer  limiter  resolves  the  discon¬ 
tinuity  slightly  better  then  with  Roe’s  fa  limiter,  but  the  discontinuity  profiles  in  both 
cases  are  sharp  (having  0(h)  width).  The  first  order  schemes  are  obviously  inferior  to 
the  second  order  schemes  in  terms  of  the  discontinuity  resolution. 
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As  follows  from  the  argument  in  Chap.2,  the  two  factors  which  determine  the  error 
in  discontinuity  location  are:  the  smooth  region  error  at  the  points  L  and  R  (see  Figure 
2.3)  and  the  error  of  extrapolation.  Sines  the  solution  in  this  example  is  a  piecewise 
constant,  points  L  and  R  can  be  chosen  far  enough  from  the  discontinuity  so  that  the 
numerical  fluxes  at  these  points  will  be  exact.  The  extrapolation  by  constant  is  exact  in 
this  case  as  well.  Therefore,  there  should  be  a  zero  error  in  the  discontinuity  location. 
Indeed,  we  have  observed,  that  the  discontinuity  location  in  the  limit-solutions  in  these 
cases  can  be  recovered  upto  the  round-off  error. 

6.2  Nonlinear  equation 

We  consider  here  the  following  version  of  Eq.(6.1): 

-  (0+)Au  +  (u2/2),  +  Uy  —  a.  (6.11) 

This  equation  can  be  rewritten  in  the  quasilinear  form 

—  (0+)Au  +  uti*  +  uy  =  a.  (6.12) 

It  is  easy  to  see  form  (6.12)  that  the  characteristic  directions  are  given  in  this  case  by  the 
vector  (ti,  1),  i.e.  it  depends  on  the  solution  (For  the  description  of  the  possible  solutions 
see  in  Sec.2.1.2). 

6.2.1  Smooth  solution 

In  case  the  boundary  conditions  are  given  by  (6.8)  and 

3  =  cos(x  +  y)(l  +  sin(®  +  y))  (6.13) 

it  is  easy  to  see  that  (6.8)  is  also  the  exact  solution  of  (6.11). 

Although  we  choose  boundary  conditions  so  that  there  will  be  no  boundary  layers 
created,  numerical  boundary  layers  may  appear.  This  is  because  some  of  the  difference 
schemes  we  shall  use  are  not  exactly  upstream.  In  order  to  avoid  any  influence  of  potential 
numerical  boundary  layers  in  our  study  of  interior  behavior,  the  error  was  measured  only 
in  the  subdomain 

O'  *  {(*,») :  1/2  <  *  <  5/2,  0  <  y  <  4/3}  (6.14) 

Numerical  experiments  with  this  model  problem  are  presented  in  Table  6.4. 

Again  the  algorithms  employing  upstream  and  homogeneous  N  schemes  produce  first 
order  accurate  solutions,  as  can  be  seen  from  Columns  1,2.  The  inhomogeneous  N 
scheme  seems  to  demonstrate  higher  order  convergence  on  coarser  grids  because  of  the 
same  reason  as  in  the  case  of  linear  inhomogeneous  problem.  The  IF  MG  based  upon 
the  Si  and  S2  scheme  leads  to  second  order  accurate  solutions  (see  Columns  4-7).  The 
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1  Difference 
|  scheme 

Upstr 

N 

horn. 

N 

inborn. 

SI 

SI 

S2 

S2 

S1&2  1 

1  Limiter 

- 

- 

- 

*!>vl 

lf>VL 

V* 

mm 

gf 

.189 

.299- 

-  .159  - 

.158 

1 .158 
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Table  6.4:  Nonlinear  inhomogeneous  problem;  slowly  varying  solution. 


reason  for  the  larger  errors  demonstrated  by  the  S2  scheme  is  the  additional  viscosity 
(the  e  correction)  which  becomes  non-zero  when  the  characteristic  directions  are  non¬ 
constant.  This  additional  viscosity  causes  that  the  S2  scheme  demonstrates  errors  on  the 
coarse  grids  even  larger  than  those  obtained  by  the  first  order  accurate  inhomogeneous 
N  scheme.  Note  that  the  choice  of  the  limiter  in  this  problem  does  not  affect  the  results. 
Column  8  presents  the  results  obtained  by  the  “blended”  S1&2  scheme  which  contains  the 
SI  scheme  with  the  weight  9/10  and  the  S2  scheme  with  the  weight  1/10.  This  scheme 
demonstrates  the  second  order  convergence  as  each  one  of  its  components  separatly, 
however,  its  errors  are  much  closer  to  those  of  the  SI  scheme  than  of  the  S2  scheme. 
Note  that  the  S1&2  scheme  when  used  with  Roe’s  Vb  limiter  is  monotonic  in  case  of  a 
homogeneous  problem. 

Table  6.5  presents  numerical  experiments  with  Eq.(6.11)  and  with  the  right-hand-side 

s  =  3  cos(3(x  +  y))(l  +  sin(3(x  +  y)))  (6.15) 

and  with  the  boundary  conditions 

u  «  sin(3(x  +  y)),  (6.16) 

which  is  the  exact  solution  of  this  problem  as  well. 

Again  the  upstream  and  homogeneous  N  scheme  demonstrate  the  first  order  con¬ 
vergence  for  this  problem  (Columns  1,2)  The  inhomogeneous  N  scheme  leads  to  much 
smaller  errors  and  seems  to  be  higher  order  accurate  on  the  coarsers  grids  (Column  3). 
The  SI  scheme  starts  to  demonstrate  second  order  convergence  on  the  finer  levels,  but 
more  than  2  cycles  on  the  finest  level  achieve  a  second  order  accurate  solution  (Columns 
4,5)  (see  also  Columns  4-6  in  Table  6.2  and  Sec2.3  in  [4]).  The  S2  scheme  does  not 
dearly  demonstrate  the  second  order  convergence  even  on  the  finest  levels  used  in  tins 
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Table  6.5:  Nonlinear  inhomogeneous  problem;  oscillatory  solution. 


experiment  (Columns  6,7)  and  leads  to  much  larger  errors  than  the  Si  scheme.  The 
S1&2  scheme  starts  to  be  second  order  convergent  on  the  finest  levels  and  leads  to  errors 
slightly  larger  than  those  of  the  Si  scheme.  The  choice  of  the  limiter  does  not  affect  the 
results. 


6.2.2  Discontinuous  solution 

Consider  again  Eq.(6.11)  with  s  =  0. 


Shock  wave 

If  we  choose  (6.10)  to  be  the  boundary  conditions  for  this  problem,  it  will  be  the  exact 
solution  for  it  as  well.  Figures  6.4(a  —  c)  present  plots  of  the  numerical  solutions  to  this 
problem  by  the  2 FMG  algorithm  employing  the  S1,S2  and  S1&2  schemes  respectively. 

The  choice  of  limiter  in  this  problem  has  no  influence  on  the  results.  The  shock 
profiles  in  case  of  the  SI  and  S1&2  schemes  are  similar,  however  the  one  created  by 
the  S2  scheme  is  less  sharp  because  of  the  larger  cross-stream  viscosity  due  to  the  e- 
correction.  All  the  profiles  have  0(h)  width. 


Rarefaction  wave 

Consider  the  same  homogeneous  equation  as  before  with  boundary  conditions 


u  =  1  -  2H(x  -  1.5), 


(6.17) 


Figure  6.4:  continued;  (c)  -  S1&2  scheme. 


Figure  6.5:  Rarefaction  wave;  (a)  -  non-physical  solution,  (b)  -  Si  scheme  with  Roe’s 
limiter. 


6.2  Nonlinear  equation 


Figure  6.5:  continued;  (c)  -  S 2  scheme,  (d)  -  S1&2  scheme 
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where  H  is  again  the  Heaviside  function.  Its  plot  is  presented  on  Figure  6.5a.  Such 
a  solution  can  be  obtained  also  by  the  upstream,  N  or  Si  schemes  (as  well  as  by  the 
Spekrejse  scheme).  This  is  because  these  schemes  demonstrate  zero  residuals  on  such 
a  grid  function.  However,  this  solution  is  non-physical,  it  violates  the  entropy  law  - 
there  are  characteristics  which  go  out  from  the  discontinuity  and  not  from  the  boundary. 
Figure  .6*56  presents  the  solution  ito  this  problem  obtained  by  2 FMG  algorithm  based 
upon  the  SI  scheme  with  tpi  limiter.  We  can  observe  that  an  (inadmissible  discontinuity 
starts  to  develop.  Figure  6.5c  presents  .the  solution  to  this  problem  obtained  by  the 
same  algorithm  based  upon  the  S2  scheme.  (The  choice  of  the  limiter  does  not  affect  the 
result  in  this  case  as  well.l  This  solution  does  not  contain  a  non-physical  discontinuity. 
This  is  due  to  the  e  correction  which  adds  some  additional  viscosity  to  the  scheme  when 
characteristic  directions  are  not  constant. 

Figure  6.5d  presents  the  solution  obtained  by  the  2 FMG  algorithm  employing  the 
blended  S1&2  scheme.  No  discontinuity  is  created  as  well  as  in  the  case  of  the  S2  scheme. 


6.2.3  Non-constant  solution  containing  discontinuities 

In  order  to  examine  the  accuracy  of  the  method  in  smooth  regions  in  the  presence  of 
discontinuities  and  to  test  the  discontinuity  locating  technique,  we  shall  construct  a 
model  problem  with  a  known  non-constant  solution  which  contains  a  shock  at  a  known 
location.  Regarding  y  as  a  time-like  direction,  consider  the  following  hyperbolic  problem 

(uJ/2).  +  uv  =  0,  (6.18) 

with  initial  conditions  given  along  the  x  axis  by 

ujy-o  =  uo  =  .5stn(xx)  +  .5.  (6.19) 

Let  us  compute  the  exact  solution  to  this  problem.  It  is  constant  along  characteristics 
and  therefore  characteristics  are  straight  lines.  Therefore,  for  every  point  (x,y)  we  can 
find  the  point  (xo,  0)  where  the  characteristic  line  which  goes  through  it  crosses  the  x-axis 
by  solving  the  implicit  equation 


x0  =  x  —  u0(x0)y.  (6.20) 

Of  course,  the  solution  is  not  unique.  But  if  the  entropy  condition  is  imposed,  or  if 
the  non- viscous  limit  of  the  viscous  equation  (6.11)  is  considered,  the  solution  becomes 
unique.  It  is  easy  to  see  that  inside  the  rectangle  (6.2)  it  will  contain  one  shock  wave 
(see  Figure  6.6)  going  along  the  line 


y  =  2x  -  2.  (6.21) 

Consider  again  the  domain  (6.2),  and  let  (6.19)  be  the  boundary  conditions  for 
Eq.(6.11). 
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r 


Figure  6.6:  Characteristic  lines  and  shock  wave  in  the  exact  solution. 
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Table  6.6:  Nonlinear  problem;  varying  solution  containing  a  shock  wave. 
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The  numerical  experiments  presented  by  Table  6.6  and  Figure  6.7  are  performed  with 
this  model  problem.  The  2 FMG  algorithm  is  used.  In  order  to  avoid  any  influence  of 
numerical  boundary  layers  and  the  shock  we  measure  the  solution  error  L\  norm  in  the 
following  subdomain 

ft"  =  {(x,y) :  1/2  <  x  <  5/2,0  <  y  <;  2x  -  3  or  max(0,2x  -.1)  <  y  <  4/$}.  (6.22) 

The  intersection  point  between  the  shock  and  a  grid  line  parallel  to  the  x  axis  is  calculated 
according  to  the  method  described  in  Chap.2.  The  segment  LR  is  taken  to  be  6  meshsizes 
long  with  its  middle  dose  to  the  assumed  shock  location.  The  extrapolation  by  a  constant 
is  used.  The  error  in  shock  location  is  measured  for  1  <  y  <  2,  where  the  discontinuity 
is  strong  enough. 

Each  pair  of  columns  in  Table  6.6  presents  one  experiment.  The  first  column  in  each 
pair  displays  the  history  of  the  L\  norm  of  the  solution  error  and  the  second  « the  L\  norm 
of  the  shock  location  error.  Again  we  can  conclude  that  a  second  order  accurate  solution 
(in  smooth  regions  and  also  in  terms  of  the  discontinuity  location)  can  be  obtained  by 
the  2FMG  algorithm  using  either  the  Si  or  the  S2  schemes.  However,  the  SI  scheme 
leads  to  smaller  errors  than  the  S2  scheme.  The  choice  of  the  limiter  does  not  affect  the 
solution.  The  S1&2  scheme  leads  to  a  solution  error  slightly  larger  than  the  Si  scheme, 
but  to  a  smaller  error  in  the  discontinuity  location.  Figures  6.7(a  —  c)  present  plots  of 
the  numerical  solution  obtained  by  the  2 FMG  algorithm  in  the  same  cases  which  are 
presented  in  Columns  2-4  of  Table  6.6. 

The  solution  plots  corresponding  to  the  SI  scheme  with  Roe’s  V’i  (Column  1)  and 
Van  Leer’s  (Column  2)  limiters  are  un distinguishable.  Therefore,  we  omit  the  first  one. 
In  spite  of  the  fact  that  the  monotonidty  of  the  solution  is  not  guaranteed  when  the 
SI  scheme  is  used  with  Van  Leer’s  limiter,  we  do  not  observe  any  oscillations  on  Figure 
6.7a.  Also  we  can  see  that  the  S2  does  not  provide  the  shock  resolution  as  sharp  as 
the  SI  scheme.  The  2 FMG  algorithm  which  uses  the  S1&2  scheme  leads  to  the  solution 
presented  by  Figure  6.7c.  The  shock  resolution  is  comparable  with  that  of  the  Si  scheme. 


Figure  0.7:  Non-constant  solution  containing  a  shock  wave;  (a)  -  SI  scheme  with  Roe’s 
0 i  limiter,  (b)  r  S2  scheme  with  Roe’s  02  limiter. 
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(c) 


Figure  6.7:  continued;  (c)  -  S1&2  scheme  with  Roe’s  rf>\  limiter. 


6.3  Choice  of  discretization 


65 


6.3  Choice  of  discretization 

Both  Si  and  S2  (as  well  as  the  blended  S1&2  scheme)  are  second  order  accurate  schemes 
when  used  with  any  of  the  limiters  considered  above. 

■  In  case  we  want  to  obtain  smaller  error?  in  smooth  regions,  sharper  shock  resolution  - 
and -better -discontinuity  location  -  the  Sf  scheme  should  be  used.  The  Sl_  scheme  is 
.also  proven  to  produce  monotonic  solution  when  used  with  Roe’s  fa  limiter  in  case  of 
homogeneous  problems.  However,  non-physical  discontinuities  can  appear  in.  the  solution 
in  case  of  a  rarefaction  wave  aligned  with  the  grid. 

The  S2  scheme  is  proven  to  be  monotonic  with  any  of  the  limiters  considered  above. 
It  does  not  admit  non-physical  discontinuities.  However,  it  leads  to  much  larger  solution 
and  discontinuity  location  errors  and  provides  shock  resolution  inferior  to  that  of  the  Si 
scheme. 

An  interesting  possibility  is  therefore  to  use  the  blended  Sl&2scheme,  i.e.  to  average 
the  SI  and  S2  schemes  with  certain  weights.  This  allows  us  to  achieve  a  more  accurate 
solution  and  better  resolution  of  discontinuities  but  still  without  admitting  the  non¬ 
physical  ones.  The  S1&2  scheme  is  monotonic  in  case  of  homogeneous  problems  when 
used  with  Roe’s  fa  limiter.  Another  alternative  can  be  to  use  the  S2  scheme  for  the  first 
multigrid  cycle  on  each  level  and  the  Si  scheme  for  the  second  one.  The  solution  error 
and  discontinuity  resolution  in  this  case  will  be  very  similar  to  those  of  the  Si  scheme. 
However,  non-physical  discontinuities  will  not  develop. 

The  choice  of  a  limiter  does  not  influence  significantly  the  resolution  of  shocks  and 
the  solution  error  for  both  Si  and  S2  schemes.  In  case  better  resolution  of  contact 
discontinuities  is  desirable,  the  compressive  limiters  (Van  Leer’s  or  Roe’s  fa)  should 
be  used.  However,  sharp  edges  in  the  solution  can  be  created  because  of  the  artificial 
compression.  When  the  SI  scheme  is  used  with  a  compressive  limiter,  the  monotonicity 
of  the  solution  is  not  guaranteed. 

We  can  conclude  that  the  choice  of  discretization  may  depend  on  which  features  of 
the  solution  are  of  more  interest  to  us. 


Chapter  7 

Discussion  and  conclusions 


7.1  Summary 

A  fast  multigrid  solver  for  a  scalar  2D  steady-state  conservation  law  is  developed. 

Two  main  problems  were  solved  in  order  to  achieve  the  goal  of  this  work: 

1.  A  new  genuinely  2D  discretization  scheme  based  on  the  compact  “9-point  box” 
stencil  was  constructed.  This  scheme  provides  a  separation  between  treatment  of  stream- 
wise  and  cross-stream  directions.  The  artificial  viscosity  is  added  in  the  streamwise 
direction  only.  High  resolution  can  be  introduced  in  the  cross-stream  direction.  As  a 
result,  numerical  solutions  produced  by  these  schemes  are  monotone  with  sharp  reso¬ 
lution  of  oblique  discontinuities  and  a  local  relaxation  process  applied  with  this  type 
of  discretization  is  stable.  The  solutions  also  demonstrate  second  order  convergence  in 
smooth  regions  in  the  case  of  a  homogeneous  problem.  This  property  is  extended  to 
inhomogeneous  equations  by  a  simple  weighting  of  the  right-hand  side. 

2.  We  have  shown  that  a  captured-shock  solution  (provided  the  discrete  scheme 
employee  is  conservative)  contains  information  about  the  discontinuity  location  up  to 
the  same  order  of  accuracy  which  is  actually  achieved  in  the  smooth  regions.  This 
assertion  has  two  important  implications.  The  first  is  the  possibility  of  performing  a 
post-pro.  easing:  once  the  discontinuity  is  recognized  in  the  solution,  its  location  can 
be  extracted  with  a  higher-order  accuracy.  The  second  implication  is  that  the  usual 
multigriJ  correction  interpolation  (bilinear  etc.)  can  be  used  in  the  neighborhood  of 
discontinuities  as  well  as  in  the  smooth  regions.  Indeed,  it  expresses  the  correct  movement 
of  the  di  continuity  provided  a  conservative  residual  transfer  (Full  Weighting)  is  used. 

Numerical  experiments  confirm  that  second  order  accurate  (both  in  smooth  regions 
and  discontinuity  location)  solutions  can  be  usually  obtained  by  the  2FMG  algorithm, 
even  with  direction-free  relaxation.  The  limit  solutions  can  also  be  obtained  by  just  few 
downstream  relaxation  sweeps  on  the  finest  level. 


7.2  Remark  on  double  discretization 

Our  original  intention  was  to  use  the  double  discretization  technique,  i.e.  to  employ  two 
different  schemes  in  the-  relaxation  and  in  the  residual  calculation.  The  scheme  used  in 
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the  relaxation  process  must  be  stable,  but  may  be  low  order  accurate.  The  scheme  used 
for  the  residual  calculation  must  be  higher  order  accurate,  but  may  be  unstable.  This 
approach  is  usually  good  for  smooth  solutions.  However,  there  is  a  certain  difficulty  to 
apply  it  in  the  presence  of  discontinuities.  Two  different  schemes  lead  usually  to  solutions 
whith  have  different  discontinuity  profiles.  Therefore,  the  residuals  calculated  by  means 
of  one  scheme  on  the 'solution  approximation  obtained  by  relaxing  another  scheme  may 
attain  large  values  in  the  neighborhood  of  discontinuities.  If  we  transfer  these  large 
residuals  together  with  others  to  the  coarser  grid,  the  numerical  process  may  become 
divergent.  However,  in  case  both  difference  schemes  used  are  conservative  the  summation 
of  the  residuals  inside  the  transition  layer  along  gridlines  crossing  the  discontinuity  will 
give  small  numbers.  Therefore,  if  this  “residual  cancelling”  is  performed  each  time  before 
going  to  the  coarser  grid,  the  multigrid  algorithm  employing  the  double  discretization 
will  produce  higher  order  accurate  solution  in  smooth  regions  as  well  as  higher  order  error 
in  discontinuity  location.  However,  such  an  algorithm  requires  to  detect  transition  layers 
representing  discontinuities  in  the  solution  and  to  perform  residual  cancelling  accross 
these  layers  on  each  grid  before  transfering  residuals  to  the  coarser  one.  This  is  the 
reason  why  we  decided  to  use  the  same  higher  order  accurate  and  stable  scheme  both 
for  the  relaxation  and  the  residual  calculation. 


7.3  Extension  to  Euler  equations 

Our  next  objective  is  to  extend  these  methods  to  the  steady-state  Euler  equations.  This 
is  not  expected  to  be  complicated,  because  a  discretization  of  the  convection  operator 
is  the  main  concern  of  the  approach  of  [3]  for  the  case  of  this  system  too.  Moreover, 
it  seems  to  be  possible  to  design  a  stable  local  relaxation  process  for  the  Euler  system, 
based  on  the  ideas  of  [3]. 


7.4  Efficiency  comparison 

We  shall  compare  now  the  number  of  operations  neccessary  to  perform  in  order  to  eval¬ 
uate  a  residual  at  one  grid  point  of  the  Spekreijse  scheme  and  of  the  Si  and  S2  schemes 
developed  in  this  work.  The  results  of  the  comparison  are  summarized  in  Table  7.1. 

We  can  conclude,  that  the  SI  and  the  Spekreijse  schemes  are  comparable  in  terms  of 
efficiency  for  residual  calculation.  The  S2  scheme  is  more  expensive.  However,  we  want 
to  note  the  following: 

•  The  schemes  constructed  in  this  work  are  based  on  narrow  stencils.  Therefore, 
they  introduce  smaller  cross-stream  viscosity  than  the  schemes  of  the  same  order 
of  accuracy  based  on  the  “dimensional  splitting”  approach  (see  Sec.2  in  [4]).  This 
means  that  the  solution  of  a  certain  quality  can  be  obtained  by  our  schemes  on  a 
grid  coarser  than  the  one  needed  for  the  Spekreijse  scheme. 
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Table  7.1:  Efficiency  comparison:  after  the  “+”  sign  -  the  number  of  additional  opera¬ 
tions  needed  to  achieve  second  order  accuracy  in  case  of  an  inhomogeneous  problem;  in 
the  parenthesis  -  the  number  of  operations  which  may  be  required  in  the  neighborhood 
of  discontinuities;  otherwise  given  is  the  number  of  operations  required  in  the  smooth 
regions. 
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•  A  pointwise  relaxation  with  the  SI  and  S2  schemes  is  stable.  To  perform  a  Newton 
iteration  at  one  point  requires  (in  addition  to  the  residual  calculation)  the  evalu¬ 
ation  of  the  numerical  fluxes  derivatives  with  respect  to  the  central  point  value. 
Tins  is  much  less  expensive  than  the  residual  calculation  (especially  if  the  approx¬ 
imate  formulas  for  derivatives  are  used).  The  evaluation  of  the  numerical  fluxes 

'derivatives  in  case  of  the  Spekreijse  scheme  is  approximately-  as  .expensive  ap  in 
case  of  the  Si  and  S2~ schemes.  However,  a  4-point  block  relaxation  is  needed  when 
relaxing  the  Spekrejse  scheme  (see.  [20]).  This  means,  that  one  Newton  iteration 
should  be  performed  . for  a  system  of  4  .  nonlinear  equations  each  time,  Instead  of 
being  applied  for  4  separate  equations  as  in  the  Si  or  S2  schemes.  This  requires  to 
invert  a  4  x  4  matrix  for  every  4  gridpoints.  This  matrix  contains  at  least  4  zero 
elements,  and  upto  8  in  case  of  a  smooth  solution.  This  additional  matrix  inversion 
increases  the  computational  cost  of  the  Spekreijse  scheme  when  comparing  to  the 
Si  and  S2  schemes. 

•  The  really  important  point  here  is  that  there  was  not  found  any  (neither  local 
nor  global)  stable  relaxation  which  can  be  applied  with  the  Spekreijse  scheme  in 
case  of  the  steady-state  Euler  system.  Therefore,  it  was  possible  to  achieve  the 
second  order  accuracy  only  by  employing  a  defect  correction  method  (see  [19]), 
which  is  not  a  fully-efiicient  way  to  deal  with  non-elliptic  problems  (see  Sec.2  in 
[4]).  Indeed,  many  cycles  were  required  for  the  Spekreijse  scheme  to  reach  second 
order  accuracy. 

We  can  conclude  that  the  multigrid  solver  for  a  scalar  2D  conservation  law  based  on 
the  discretization  schemes  developed  here  is  slightly  more  efficient  than  the  one  based  on 
the  Spekreijse  scheme.  However,  when  extended  to  the  Euler  system  they  are  expected 
to  be  much  more  efficient  than  the  solver  developed  in  [19]. 


7,5  Some  properties  and  future  development 

The  versions  of  the  Si  and  S2  schemes  which  were  shown  to  be  monotonic  for  the  case 
of  a  homogeneous  equation  seem  to  be  TVD  (total  variation  diminishing).  Although 
we  do  not  prove  this,  it  must  follow  from  the  monotonic  property  for  these  schemes 
exactly  as  in  the  ID  case.  This  means  that,  when  these  schemes  are  used  to  solve  a 
time-dependent  problem,  convergence  of  the  solution  is  theoretically  justified  for  a  finite 
time  interval.  This  does  not  contradict  the  result  of  [9],  because  when  used  to  discretize 
a  time-dependent  problem,  both  the  Si  and  the  S2  schemes  will  retain  only  first  order 
spatial  accuracy.  However,  discontinuities  will  be  sharply  resolved  in  the  solution.  The 
second  order  accuracy  can  possibly  be  obtained  by  treating  the  time  derivative  as  the 
right-hand-side  of  an  inhomogeneous  problem  (see  Sec.5),  but  such  a  discretization  will 
not  have  the  TVD  property  anymore. 

Note  that  the  use  of  TVD  schemes  has  a  theoretical  advantage  for  time-dependent 
problems  and  finite  time  evolution  only.  However,  when  constructing  a  certain  scheme 
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for  practical  purposes  it  may  be  enough  to  require  that  it  can  be  rewritten  as  to  another 
scheme  which  is  not  necessarily  of  the  positive  type,  i.e.  it  can  have  small  negative 
coefficients  (0(h)  compared  with  the  positive  coefficients)  at  some  points. 

We  can  also  distinguish  between  smooth  regions  and  the  neighborhood  of  disconti¬ 
nuities,  and  use  different  schemes  for  each  of  them:  a  more  complicated  scheme  with 
limiters  near  discontinuities  only,  and  i  simplified  scheme  in  the  smooth  regions.  These 
simplifications  can  lead  to  a  substantial  computational  cosir  reduction  of  the  method  in 
smooth  regions.  Also,  additional  local  relaxation  sweeps  may  be  performed  in  the  neigh¬ 
borhood  of  strong  discontinuities,  increasing  the  efficiency  of  the  method  for  little  extra 
cost. 

It  must  be  possible  to  construct  a  third-order  accurate  upstream  scheme  based  on 
the  same  “9- point  box”  stencil:  4  upstream  points  are  enough  to  obtain  an  0(h3)  cross¬ 
stream  truncation  error,  while  the  third  order  accuracy  in  other  directions  can  be  achieved 
again  by  a  certain  weighting  of  the  right-hand-side  of  the  equation.  An  approach  sim¬ 
ilar  to  ENO  (essentually  non-osdllatory)  in  [13]  can  be  introduced  in  the  cross-stream 
direction.  However,  some  technical  details  still  have  to  be  verified. 

Another  interesting  possibility  is  to  construct  upstream  difference  schemes  based 
on  thiner,  but  longer  stencils.  This  approach  will  improve  resolution  of  characteristic 
components  and  discontinuities  and  can  possibly  lead  to  a  very  high  order  accuracy 
compact  schemes. 
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